Homogenization of Maxwell's Equations in Periodic Composites: Boundary Effects 

and Dispersion Relations 



Vadim A. Markel * 

Department of Radiology and Graduate Group in Applied Mathematics and Computational Science, 
University of Pennsylvania, Philadelphia, PA 19104 

John C. Schotland t 

Department of Mathematics, University of Michigan, Ann Arbor, MI 48109 

(Dated: May 24, 2012) 

We consider the problem of homogenizing the Maxwell equations for periodic composites. The 
analysis is based on Bloch-Floquet theory. We calculate explicitly the reflection coefficient for a 
half-space, and derive and implement a computationally-efficient continued-fraction expansion for 
the effective permittivity. Our results are illustrated by numerical computations for the case of 
two-dimensional systems. The homogenization theory of this paper is designed to predict various 
physically-measurable quantities rather than to simply approximate certain coefficients in a PDE. 
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I. INTRODUCTION 



Theories of electromagnetic homogenization of com- 
posite materials — also known as effective medium theo- 
ries (EMTs) — have a history which dates to the time of 
J.C. Maxwell. Nevertheless, these theories continue to 
attract attention and even controversy, as evidenced by 
recent reviews [2[ and many references therein. In ap- 
plied mathematics, the theory of homogenization based 
on multiscale analysis of partial differential equations is 
also well-established [3-6]. However, interest in EMTs 
has been steadily on the rise for the past ten years with 
conceptually new approaches continuing to appear 0- 
[Tq| . This can be explained, perhaps, by noting that 
the tasks of relating the existing mathematical theories 
to physical observables and of determining the range of 
applicability of a given theory have not been fully ad- 
dressed, particularly, for the case of Maxwell's equations. 
Indeed, in the past ten years or so, homogenization theo- 
ries have been applied to obtain "extreme" properties of 
electromagnetic composites, including the phenomenon 
of strong "artificial" magnetism. At the same time, a sig- 
nificant experimental progress has been recently made in 
manufacturing deeply- subwavelength (in the visible spec- 
tral range) periodic metallic nanostructures [lll-[l3| . The 
question is whether the existing theories are directly ap- 
plicable or accurate enough to guide the experimental 
design of periodic nanostructures of desirable properties. 
Another reason for the renewed interest in homogeniza- 
tion theories is that, in addition to abstract mathemati- 
cal results, there is a need for efficient, stable computa- 
tional methods. Thus the question of how to construct 
physically-relevant and computationally-effective EMTs 
and determine their limits of applicability have not been 
completely settled. 
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This paper is an attempt to address the above issues 
for the case of periodic composites; random media are 
not considered. The framework we develop is based on 
the Bloch-Floquet expansion, which is a well-known tool 
in homogenization theory [5l ll4Ml9| . However, in several 
aspects, we go beyond the standard theory. In partic- 
ular, (i) we explicitly account for boundary effects and 
derive a general expression for the reflection coefficient 
(many existing homogenization theories consider infinite 
composites) (ii) we make use of the integral equation 
formulation of scattering theory for the Maxwell equa- 
tions. The resulting formulas for the effective medium 
parameters (EMPs) have a different mathematical struc- 
ture than those derived from partial differential equa- 
tions (iii) we develop a computationally-efficient algo- 
rithm for calculating the EMPs. The algorithm is based 
on a continued-fraction expansion of the self-energy and 
is obtained from a new result on the resolvent of a lin- 
ear operator and (iv) a numerical study of stability and 
convergence is performed for some test cases. Stability 
is investigated by comparing the results for inclusions of 
the same volume fraction but different shape and of the 
same shape but different volume fractions. 

It is useful to recognize that all EMTs can be clas- 
sified as either standard or extended. A standard EMT 
is obtained by taking the limit h 0, where h is the 
scale of the medium's heterogeneity; in this paper, h is 
the lattice spacing. In standard theories, h is viewed as a 
mathematically- and physically-independent variable and 
the resulting EMPs are independent of h, as long as the 
latter is small enough for the theory to be applicable. An- 
other feature of all standard theories is the so-called law 
of unaltered ratios [20|, which states that, if a composite 
medium is made of several constituents with permittiv- 
ities €j (j = 1,2,...) and if ej — >• Ae^ (A > 0), then the 
effective permittivity e also scales as Ae. 

Extended EMTs came to the fore (at least in the 
physics literature) in [2l|, [22|. The basic idea of these 
papers is to note that one can compute the exact electric 
and magnetic polarizabilities, a e and a m , of a spherical 
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particle through the use of the first Lorenz-Mie coeffi- 
cients, a\ and &i, even when the sphere in question is 
not small compared to the external wavelength. These 
polarizabilities can be used to construct an "extended" 
Maxwell- Garnett approximation. Since a\ and b\ are not 
proportional to the sphere volume, except in the qua- 
sistatic limit, the resultant EMTs contain the sphere ra- 
dius explicitly. In Refs. [23|, |24[, a counter-intuitive effect 
of non-commuting limits was demonstrated. Specifically, 
it was shown that insofar as the effective refractive index 
of a photonic crystal is computed from the slope of the 
dispersion curve near the T-point, different results are 
generally obtained depending on which of the two limits, 
h —> and e\ — > oo is taken first, where e\ is the permit- 
tivity of one of the components of the photonic crystal. A 
related point is that a complete theory of homogenization 
requires error estimates. That is, it is essential to deter- 
mine how the error in the homogenization limit depends 
upon contrast. Moreover, the reflection and transmis- 
sion properties of the composite medium have not been 
considered [HHII. 

In this paper, we develop a standard EMT. However, 
when considering reflection and refraction at a planar 
interface, we derive formulas for the reflection and trans- 
mission coefficients, which are valid for finite values of 
h. Then we show that taking the limit h results in 
the standard Fresnel coefficients. In this case, the elec- 
tric and magnetic properties of the medium constituents 
do not mix, in agreement with [25[. That is, if we begin 
with nonmagnetic inclusions, the resultant composite is 
also nonmagnetic. An extended EMT can be obtained 
by taking a different limit, in which the permittivity of 
one of the constituents scales as 1/h 2 [26[ . Here we note 
again the existence of the effect of non-commuting lim- 
its [23|, [24I H3 -> which calls for additional scrutiny of the 
homogenization results thus obtained. In particular, one 
would expect that, in the limit considered in [26|, Fresnel 
formulas would also be reproduced, but with a nontriv- 
ial magnetic permeability. We have not been able to 
show that this is the case. In other words, it is not clear 
whether the EMPs obtained from an extended EMT are 
independent of the incidence angle or, more generall y, o f 
the type of incident wave. This is in accord with [28|-[33|, 
which find that the conditions under which metamate- 
rials exhibiting strong magnetic resonances can be as- 
signed purely local (incidence-angle-independent) EMPs 
are rather restrictive. The same point has been made in 
the recent review article [2|. 

An additional feature by which EMTs can be classi- 
fied is the physical model of the medium. In the model 
of dipole lattices, the medium is thought of as being com- 
posed of point particles which are completely character- 
ized by their polarizabilities (electric and, possibly, mag- 
netic) and whose shape and size do not enter into the 
problem directly [3J-[36| . Alternatively, one can consider 
the space as a two-component continuous medium [371 - 
[39|. The point-dipole model is appealing because of its 
simplicity but leads to serious mathematical problems. 
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FIG. 1: (color online) Sketch of the geometry considered: an 
infinite 3D lattice. 

The so-called dipole sum (also known as the lattice sum 
or the dipole self-energy), which plays a key role in this 
model, diverges in the case of three-dimensional lattices. 
While it is true that even divergent series can be summed 
by means of applying various mathematical tricks, the 
results obtained depend on the particular trick used, a 
state of affairs that is not very satisfying. Therefore, we 
will adopt from the start a model of a two-component 
continuous medium. As the development in this paper 
progresses, it will become apparent why the point-dipole 
model is inadequate. 

The mathematical development in this paper begins 
by considering the integral equation obeyed by the po- 
larization field, which is introduced in Sec. [Ill In Sec. IIII( 
we derive a homogenization theory of the standard type 
for infinite periodic media. Reflection and refraction at 
a planar boundary is considered in Sec. IIVI In Sec. |VJ 
we discuss the correspondence between the point-dipole 
model and the continuous-medium model of this paper. 
One mathematically-novel element of the theory devel- 
oped herein is a continued-fraction expansion of the ef- 
fective permittivity, which is derived in Sec. IVII and used 
in the numerical simulations of Sec. IVIII The expansion 
has its origins in a theorem on resolvents of general linear 
operators (with no special symmetry properties), which 
is stated in Sec. IVII and proved in the appendices. A dis- 
cussion and a brief summary of results are contained in 
Sees. IVHU and ED 



II. BASIC EQUATIONS 

The geometry of the problem we consider is sketched 
in Fig. [TJ The medium consists of two intrinsically non- 
magnetic constituents: a host medium of permittivity 
and periodically- arranged inclusions of permittivity e a . 
In practice, the host is often a transparent dielectric with 
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Ree^ > 0, < Ime^ <C Ree^, and the inclusions are metal- 
lic. However, the theory of this paper places no such 
restriction on the permittivities and only requires that 
Ime^ > 0, Ime a > 0. In the case when the host medium 
is vacuum, we will take e& = 1 + i0. The inclusions are 
arranged on a cubic lattice of period h. The position vec- 
tor of the center of each unit cell is denoted by r n , where 
n can be viewed as a composite index: n = (n x ,n y ,n z ) 
and r n = h(xn x + yn y + zn z ). Whenever a summation 
over n (or a similar composite index m) appears in the 
text, it is implied that the sum runs over all three in- 
teger indexes. Inside the nth cell, the spatial region Q n 
has the permittivity e a , and the rest of the cell has the 
background permittivity e&. All regions Q n are identical 
and only differ by translation. It is assumed that Q n can 
touch but not cross the cell boundaries. No assumption 
on the connectivity of Q n is made. The union of all re- 
gions Q n is denoted by fi to t an d the volume of each region 
by V: 



and 



fit. 



IJfi" - 



d r = V . 



(1) 



We work in the frequency domain and the common factor 
exp(— iujt) is suppressed. All frequency-dependent quan- 
tities, such as the permittivities and e a , are evaluated 
at the frequency uj. 

The mathematical development in this paper begins 
with the integral equation 



Ei(r)+ f G(r,r')P(r')dV 



, r e fit 



Here P(r) is the vector of "polarization" , which is related 
to the electric field E(r) by 



P(r, - 1^E W 



(3) 



Ei(r) is the incident electric field, G(r, r') is the regular 
part of the free-space, retarded Green's tensor, and 



X 



2e 6 



(4) 



Note that P(r) defined in (j3j) is not the true physical po- 
larization, which is given by [e(r) — l]E(r)/47r, but rather 
an auxiliary field; P(r) vanishes in the host medium while 
the true polarization does not. 

In what follows, we will make use of the spatial Fourier 
transform of the Green's tensor, namely, 



G(r,r') 



4-7T 



/ 



d 3 p 



if(p)exp[ip-(r-r')] , (5) 



where 



K(p) = 



2k 2 b +p 2 



3p i 



V 2 ~ k 2 



(6) 



k^j — ejjk , k 



c 



(7) 



Here the wave number in the background medium is de- 
noted by kb and the wave number in vacuum by k. We 
note that the integral equation (j2]) is equivalent to the 
pair of curl Maxwell equations written in the frequency 
domain. 



III. WAVES IN INFINITE LATTICES 

A. Three-dimensional lattices 

Consider the propagation of a wave in a three- 
dimensional infinite lattice. In this case, the incident 
field is absent and Eq. (J2j) must be satisfied for E^ = 0. 
We seek the solution to Eq. (|2]) in the form of a Bloch 
wave: 

P(r) =exp(iq-r n )F(r-r n ) , r e Q n . (8) 

Here q is the Bloch wave number and F(r) is a vector 
function. Equivalently, if we write r = r n + R, then 

P(r n + R) = exp (iq • r n ) F(R) , RgD. (9) 

In this formula, Q = Qo is the region centered at the 
origin of a rectangular reference frame. From the above 
relation, we find the equation obeyed by F(R): 

F(R) = ^ / WCR^FiTL'^R' , (10) 

where 

W(R, R') = G ( r n + R, r m + B!) 

m 

x exp [iq • (r m - r n )] . (11) 

It can be seen that W is independent of n. It should also 
be noted that the summation in Eq. ([TT]) runs over the 
entire lattice, including the term m = n. In theories that 
consider point-like particles, the dipole sum is defined 
as an incomplete lattice sum, which excludes the term 
m = n. This makes application of the Poisson summation 
formula problematic and unnecessarily complicates the 
mathematics [36|. 

Returning to our derivation, we evaluate W as 



w(R, R') = y / ?S^(p) ex P fa> • ( 



(R - R ) 



4ir 
3/? 



(2tt) 3 

x ^2 exp [i(p - q) • (r„ - r m )] 

m 

^X(q + g)exp[i(q + g).(R-R')] , (12) 



4 



where 



2tt , a 

— (xn, + y% 



■zn z ) 



(13) 



are the reciprocal lattice vectors and we have used the 
Poisson summation formula 

^exp[i(p-q) • (r m - r n )] = (^-j X^( p_q_g ) " 

(14) 

The summation in Eqs. ([T2]) , ([H]) is over the complete set 
of reciprocal lattice vectors; equivalently, it can be viewed 
as summation over the triplet of indexes (n X: n y ,n z ) 
which appear in ([T3]) . 

The series in the right-hand side of ([T2]) diverges when 
R = R/. This is the well-known divergence of the dipole 
sum [40j which hinders the analysis of waves in lattices 
made of point-like polarizable particles. The model of 
point-like dipoles is discussed in more detail in Sec. [Vj In 
the equations derived above, the divergence is of no con- 
cern because W(R, R') appears only inside an integral 
and the singularity in question is integrable. 

Upon substitution of ([T2]) into (fT0]h we obtain 

F ( R ) ^E^+S) ex P [*(<! + S) • R ] 
g 

x / F(R / )exp[-z(q + g) • R'] • (15) 
It follows from ([T5]) that F(R) can be expanded as 

F(R) = F g exp [i(q + g) • R] (16) 



and that the expansion coefficients satisfy the system of 
equations 



(17) 



where p = V/h 3 is the volume fraction of inclusions and 
M(g) is defined by the expression 



Af(g) = ^ / exp(-ig-R)d 3 i? 



(18) 



Note that M(g) is defined only by the shape of the in- 
clusions and is invariant with respect to the coordinate 
rescaling r — ^ Ar. Some mathematical properties and 
calculations of M(g) for special geometries are given in 
Appendix [A] 

So far, we have simply restated the well known theo- 
rem of Floquet. The eigenproblem ([T7]) defines the band 
structure of a photonic crystal. It is well known that 
EMTs are not always applicable to photonic crystals. 
However, there exists a regime in which EMPs can be 
reasonably introduced, and this regime will be explored 
below. Namely, if qh, k^h <C 1, we can consider the cases 



g = and g ^ in (fT7|) separately. This yields the 
following equations: 



F = p X K(q) 



F g = PxQ(g) 



g)F g 



M(g)F +^M(g-g')F g 

g'^0 



(19a) 
(19b) 



where 



Q(s) 



lim K(q + g) = 1 - 3g ® g , g + . (20) 

h— >0 



Here g = g/|g| is a unit vector. 

The derivation of Eqs. (fT9j) is one of the key develop- 
ments of this paper. It can be seen that the equations in 
(|19b|) do not contain the variables k or q, but are com- 
pletely defined by the geometry of inclusions and by the 
variable %. Moreover, these equations are invariant with 
respect to the rescaling r — » Ar. For any given value of 
Fo, (|19bp can be solved uniquely as F g = A g Fo, where 
the tensors A s depend on g, the shape of inclusions, and 
on X- Given this result, we can write 



-g)A g F = SF 



(21) 



where the tensor E has all the properties of A s and, in 
addition, is independent of g. It will be shown in Sec. IVII 
that E plays the role of the self-energy and originates 
due to the electromagnetic interaction within and be- 
tween the inclusions. It will also be shown that E can 
be computed as a resolvent of a linear operator, which 
depends only on the shape of inclusions. 

Using the notation introduced in (f2T]h we can rewrite 
Jal) as 



(22) 



[i-^(q)(i + s)] f = o . 

This equation has nontrivial solutions if 

det [l-p X Jf(q)(l + E)] = 0. 



(23) 



Here the quantity in the square brackets is a 3 x 3 matrix. 
For a fixed value of k (that is, at a fixed frequency), the 
condition ([23]) is an algebraic equation with respect to 
the Cartesian components of the Bloch vector q. Roots 
of this equation, computed at different values of /c, de- 
termine the dispersion relation q(fc). There can be more 
than one branch of the dispersion relation corresponding 
to different polarization states. By polarization of the 
mode, we mean here the direction of the vector Fo. 

EMPs can be inferred by comparing these results to 
the polarization states and dispersion relation in a ho- 
mogeneous medium characterized by tensor permittivity 
and permeability e and ft. However, it is not possible 
to determine e and ft simultaneously and uniquely from 
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consideration of the dispersion relation alone. For exam- 
ple, in an isotropic medium, only the product of these 
two quantities (the squared refractive index) can be un- 
ambiguously obtained. Indeed, the dispersion relation in 
such a medium is invariant with respect to the transfor- 
mation e — > £e, ft — » where £ ^ is a complex 
number. To determine e and Jl uniquely, one must con- 
sider reflection and refraction at the medium boundary. 
This will be done in Sec. [iVl In particular, it will be 
shown that, in order to obtain the correct Fresnel reflec- 
tion coefficients, one must set Jl = 1. 

To summarize the results of this section, the electro- 
magnetic modes of a medium can be found if the tensor E 
is known. Computation of the modes involves diagonal- 
ization of a 3 x 3 matrix, while the tensor E is uniquely 
determined by the solution to Eqs. (|19b|) . The latter is 
an infinite set of equations which must be appropriately 
truncated in numerical computations. Thus, we have re- 
duced the homogenization problem to solving a set of 
algebraic equations in which the shape of the inclusions 
appears only in the functions M(g). 



B. Main homogenization result for 
three-dimensional composites with well-defined 
optical axes 

The standard description of electromagnetic waves in 
anisotropic crystals is based on the assumption that the 
tensors e and jl commute and are simultaneously diago- 
nalizable by a rotation of the reference frame, with purely 
real Euler angles. The axes of the reference frame in 
which e and Jl are diagonal are known as the optical 
axes. Moreover, standard textbooks often specialize to 
the case Jl = 1 which is a very good approximation in 
crystal optics [41]. In the most general case, however, 
the tensors e and Jl do not commute, which gives rise 
to two distinct sets of electric and magnetic axes. Fur- 
thermore, e and Jl are complex- valued, symmetric and, 
hence, non-Hermitian matrices. A purely real rotation 
that diagonalizes any one of these two tensors may not 
exist. A mathematically tractable dispersion relation for 
the most general case has been derived only recently [42| , 
and we will use below one particular case of this result. 

For the composite medium consisting of non-magnetic 
components, which is considered in this paper, the situ- 
ation is somewhat simpler. It can be seen from Eq. ([23]) 
that a unique set of optical axes exists if the tensor E 
is diagonalizable by a real-angle rotation of the reference 
frame. Thus, the issue of commut ability of two different 
tensors does not arise in this case. 

In this subsection, we assume that the optical axes of 
the composite medium (that is, the principal axes of the 
tensor E) exist and, moreover, coincide with the crys- 
tallographic axes of the medium. The latter assumption 
is not really necessary but any composite can be cut is 
such a way that it holds. In particular, E is diagonal in 
the reference frame defined by the crystallographic axes 



(which is the laboratory frame in this paper) if the in- 
clusions are symmetric with respect to reflections in each 
of the xy-, xz- and yz-planes. The principal values of E, 
denoted by S aa (a = x,y,z), are not necessarily equal 
in this case. The two familiar examples of reflection- 
symmetric inclusions which result in all three principal 
values being different are a general parallelepiped and an 
ellipsoid with unequal semi-axes. However, if the inclu- 
sions also have cubic symmetry (which, in addition to 
reflections, includes rotations about each axis by the an- 
gle 7r/4), then E is reduced to a scalar and the effective 
medium is isotropic. 



1. General direction of propagation 

Let the tensor E be diagonal in the rectangular frame 
xyz. We then use the expression (j6]) for if(q), evaluate 
the determinant in Eq. (j23j), and obtain the following 
equation: 



4Ug [1 + 2^(1 + Sea)] , 



= (fc,q) = 0, (24) 



where 



and 



: (fc,q) = k 4 - £/ c {q)k 2 + <% c (q) 



(25) 



^c(q) =ql ( — + — )+ Q 2 y (— + ~ 

\Vv Vz) y \Vx Vz 



+Q Z 



1 1 

Vx Vv 



(26a) 



*c(q) 





+■ qt + 


Q 4 Z , / 




VyVz 


VxVz 


VxVy Vz \ 


,Vx ^ 






) + &z ( 1 


1 


Vy 


\rjx Vz 


/ Vx \Vy 


+ — 

Vz 



The quantities rj a are given by 

1 + 2p X (l + Z aa ) 
1 - p X (l + E aa ) 



Va Q>- 



. (26b) 



, a = x,y,z (27) 



and the subscript in srf c and £$ c has been used to 
emphasize that these expressions are applicable to com- 
posite media and have been obtained by evaluating the 
left-hand side of (123]) . 

The set of dispersion relations (|24|) - (|26|) should be com- 
pared to the analogous set of equations in a homogeneous 
medium characterized by the effective tensors e and p. 
Generally, the dispersion relation in such media reads 

det [(q x /i _1 qx) + k 2 e] = , (28a) 

if exists, or 

det [(q x e _1 qx) + k 2 Jl] = , (28b) 
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if e _1 exists. If both jl and e are invert ible, the two 
equations (j28a|) and (|28b[) are identical. 

For homogenization theory to be applicable, the effec- 
tive medium must have the same symmetry as the com- 
posite. It is evident, therefore, that the principal axes of 
E should coincide with the optical axes of the effective 
medium. Denote the principal values of e and ft by e aa 
and ft a a- Let us further assume that ft is invertible. In 
this case, Eq. (|28a|) takes the following form: 



where 



and 



0fc(*,q) 



c£yy£zz@h(k,q) = 



4(q)fe 2 + ^(q) 



(29) 



(30) 



*4(q) =<£ 



l 



l 



tyyl^zz tzzfAyy 



l 



l 



^xxl^yy ^yyl^xx 
^4 



€xxfJ>zz ^zzf^xx 

(31a) 



#fc(q) 



^yy^zzl^yyl^zz 


^xx^zzf^xx, 


ihl ( i 




€-zzl~l>zz \^xxf^yy 




<&<& ( i 




tyyl^yy \^xx^zz 




q 2 y Qz ( i 




^xxf^xx \^yyf^zz 





(31b) 



Here the subscript in ^4 and has been used to 
emphasize that these expressions are applicable to homo- 
geneous media. In the case Jjl xx = j2 yy = p zz = 1, ([29]) 
reduces to the well-known Fresnel equation. 

The prefactors in Eqs. and ([29]) are "almost 

always" nonzero, except in the case of non-dissipative 
plasmas, which can support longitudinal waves. This 
case will be considered by us separately. Assuming that 
the prefactors are nonzero, the dispersion relations are 
@ c (k, q) = f° r the composite medium and f^(fc, q) = 
for the homogeneous medium. We can introduce EMPs 
for the composite by observing that these two dispersion 
relations become identical if we set 



(32) 



where £ ^ is an arbitrary complex number. As was 
already mentioned, the non-uniqueness in the above def- 
inition of the EMPs can not be removed by considering 
the dispersion relations alone. 

Several remarks regarding the dispersion relations ob- 
tained above should be made. First, in the general case, 
the functions @ c (k, q) and ^(fc, q) can not be factorized 
into products of two quadratic forms in the variables /c, 
q x , q y and q z . However, such a factorization becomes 
possible for special directions of propagation, when one 



or more of the Cartesian components of q are zero. Ex- 
amples will be given below. 

Second, the condition ([32]) . which guarantees that 
^ c (fc,q) = ^(fc,q), requires that the effective perme- 
ability p be a scalar. Any deviation of ft from a scalar 
will result in different laws of dispersion in the composite 
and in the effective medium with no hope of obtaining the 
same measurables from these two models. This require- 
ment that ft be a scalar even in a strongly anisotropic 
composite is difficult to justify on physical grounds, un- 
less, of course, fx = 1. 

Third, the dispersion relations ([23]) (for a composite 
medium) and (|28a|) . (]28b|) (for a homogeneous medium) 
appear to have very different mathematical structure. 
The fact that they reduce to the same equation under 
the simple condition ([32]) is quite remarkable. 

Thus, we have shown that, if orthogonal optical axes 
of the composite medium can be defined, its dispersion 
relation q(fc) and its isofrequency surfaces [defined as the 
sets containing all q such that @ c (k, q) = for each k = 
uj/c] are equivalent to those obtained in a homogeneous 
medium with EMPs e and fi given by ([32]) . where the 
quantities rj a are defined in ([27]) . 

Since it will be proved below that the correct choice 
of the parameter £ in ([32]) is £ = 1, we now state the 
main homogenization result of this paper pertaining to 
the principal values of the EMPs: 



l + 2px(l + £aa) 
l-px(l + Saa) 



Mo 



= 1 . 



(33) 



It can be seen that the Maxwell- Garnett mixing formula 
is obtained from ([33]) by setting E = 0. Electromag- 
netic interactions of inclusions in the composite result 
in a nonzero value of E and, correspondingly, in the de- 
viation of the EMPs from the predications of Maxwell- 
Garnett theory. 



2. Propagation along crystallographic axes 

Consider a plane wave propagating along the z-axis, so 
that q x = q y = 0. In this case, 



2 (fc,q) =k 4 



= k 



kVz(- + - 



Vx 



Vx 



Vy / VxVy 
Vy 



(34) 



Thus, ^ c (fc,q) is factorized into a product of two 
quadratic forms, giving rise to two branches of the disper- 
sion relation: q 2 = r] x k 2 and q% = f] y k 2 . Obviously, these 
two branches correspond to x- and ^/-polarized modes. 
It can be seen that, in agreement with ([32]) . the quanti- 
ties tJq, give the effective squared refractive index for the 
transverse modes of the composite. 

In addition to the two transverse modes, a longitu- 
dinally-polarized mode can also exist under certain condi- 
tions. A mode with an arbitrary wave number q a , which 
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propagates and is polarized along the same axis a, exists 
if and only if 



1 + 2p X (l + = . 



(35) 



Under this condition, the equality ([24]) holds, even if 
^c(£,q)^0. 

Let us consider briefly the physical conditions for exis- 
tence of the longitudinal waves. From the property (|A4|) 



(given in Appendix |A]) , it follows that lim p ^i S aa = 0. 
Consequently, the longitudinal waves exist in the high- 
density limit if 1 + 2% = 0, which is only possible if 
e a = 0. This is the well-known condition for longitudinal 
waves in non-dissipative plasma. The low-density limit 
can not be considered so easily because S aa does not ap- 
proach zero when p — > (see the Sec. IIII C|) and can, in 
fact, diverge for certain values of \- However, we can use 
the reciprocity substitution p 1 — p, e a 65 to see 
that, in the low-density limit, the condition for existence 
of the longitudinal waves is 65 = 0. Quite analogously, 
longitudinal waves can be obtained by considering the 
dispersion relation ([29]) and setting one of the principal 



3. Propagation in a crystallographic plane 

We now discuss the case when q lies in the xz-plane. 
Problems of this type can arise when one considers re- 
flection and refraction at the interface z = 0, where the 
xz-plane is the plane of incidence, as is shown in Fig. [2j 
Under the condition q y = 0, we have 



@ c (k,q) = k 4 -k 2 



1 



VyVz VxVy 



'e - ( q l - 






\Vy 


VyJ\ 





1 

Vz 

Vy 
2 



1 



1 

Vy 



(36) 



Thus, f^ c (fc,q) is factorized into a product of two 
quadratic forms, which correspond to the s- and p- 
polarized modes. 

By equating the first factor in ([36]) to zero, we obtain 
the dispersion relation for the s-polarized wave: 



Vy Vy 



(37) 



The vector Fo of the s-polarized wave is aligned with 
the 7/-axis and is, therefore, perpendicular to the plane 
of incidence. 

By equating the second factor in ([36]) to zero, we obtain 
the dispersion relation for the p-polarized wave: 



Q 2 X 



= k 2 



(38) 



We can now find the vector Fo for the p-polarized wave by 
considering the nontrivial solutions to ([22]) . It can be eas- 
ily seen that Fo lies in this case in the plane of incidence 
(its projection onto the ?/-axis is zero), and the x and 
z components of Fo, Fq x and Fq Zi satisfy the following 
relation (details of derivation are given in Appendix IB]) : 



Fox 1 + 2px(l + Vzz) Qz 
F 0z 1 + 2px(l + & 



(39) 



Eq. ([39]) will be used below in Sec. IIVI to compute the 
half-space reflection coefficient for the p-polarized inci- 
dent wave. 



C. Low-density and low-contrast limits 



Iteration of Eq. ([19b[) results in the following expansion 
for the self-energy: 

E = p X E M(-g)Q(g)M(g) + (p X ) 2 

x M (-g)Q(g) M (g-g / )Q(g / )^(g / ) + --- (40) 

g,g'^o 

It is important to note that this expansion should be 
used with caution. Indeed, if x is of the order of unity 
or larger, the series in ([40]) does not converge, even for 
arbitrarily small values of the density p. This result may 
seem unexpected, but it is easily understood by observing 
that the functions M(g) depend on p and obey the sum 
rules ([A2]) . 

In Sec. [VH a more useful (and always convergent) ex- 
pansion of E will be derived. Here we note that the 
functions M(g) are independent of x- Therefore, ([40]) is 
the formal expansion of E into the powers of x- Thus, 
in the low-contrast limit (x —> 0), we have E — ^ PX^ii 

where a\ = Zlg^o ^( _ g) ( 9(g)^(g)- In tne case °f 
three-dimensional inclusions with cubic symmetry, <j\ is 
identically zero. Then the first non- vanishing term in the 
low-contrast expansion of E is given by (px) 2(J 2, where <r 2 
grows naturally out of the second term in the right-hand 
side of (1401). 



D. Two-dimensional lattices 

Consider a medium in which e = e(x, y) is independent 
of z. As above, we assume that e(x, y) is periodic on a 
square lattice with lattice step h. The homogenization 
theory for this medium can be obtained either by consid- 
ering a three-dimensional lattice with unequal steps h Xl 
h y , h z and taking the limit h z — >• 0, or by following the 
derivations of Sec. IIII Al taking account of the modified 
geometry. The results obtained are very similar to those 
in the 3D case, with some obvious modifications. Specif- 
ically, we arrive at Eqs. (]19a|) , (|19b|) in which, however, 
we must take g = (27r//i)(xn £C + yn y ). Additionally, 
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in the integrals (fT8]h ft must be understood as a two- 
dimensional region (the intersection of an inclusion with 
the x?/-plane), V as the area of £7, and d 3 R is replaced 
by d 2 R. The definition of Q(g) J2QD remains unchanged, 
but Q(g) is now a 2 x 2 tensor. 

Consider a wave propagating in the x?/-plane and po- 
larized along the z-axis. In this case, F g = zF g , where F g 
is a scalar and E can be found analytically in general. In- 
deed, we have in this case Q(g)F g / = F g /, Q(g)Fo = Fo, 
and Eq. ([19b)) becomes 



PX 



M(g)F + ^ M(g - g')F g 

g'#0 



The solution to this equation is 

PX 



(1 - p)x 



M(g)F 



(41) 
(42) 



where some of the properties (|A2jl have been used (keep- 
ing in mind that the term g = must be excluded from 
the summation). We then have 



PX 



(i - p)x 



^M(-g)M(g) 



(1 - P)X 



- (1 - P)X 

(43) 

It can be seen from the above equation that Yt zz does not 
approach zero when p — » 0, as was discussed in Sec. lIII Cl 
Upon substitution of ([43]) into ([33]) . we find that 



(1 - p)e h + pe a = (e) . 



(44) 



Thus, the effective permittivity for z polarization is given 
by the arithmetic average of e(x, y). This is in agreement 
with Krokhin al Hill. 



E. Concept of the smooth field 

The result ([44]) for a z-polarized wave could have been 
anticipated. To understand better why the effective per- 
mittivity in this case is given by an arithmetic average, it 
is instructive to consider the concept of the smooth field. 
The smooth field S(r) changes slowly on the characteris- 
tic scale defined by the heterogeneities in the medium. As 
a result, one can factorize spatial averages of S(r) multi- 
plied by any rapidly- varying function. For example, we 
can write (Se) = (S)(e), etc. 

Let us recall some well-known results for ID 
periodically-layered media [43]. The effective permittiv- 
ity of such media is ey = (e) for waves polarized parallel 
to the layers and e± = (e -1 ) -1 for waves polarized per- 
pendicularly to the layers. These two results can be ob- 
tained quite expeditiously by applying the concept of the 
smooth field. In the case of tangential polarization, the 
electric field E is smooth. This follows from the bound- 
ary condition which requires that the tangential compo- 
nents of the electric field be continuous at all interfaces. 



Consequently, we can write 



(D) = (eE) = (e)(E) , 



(45) 



from which it follows that ey = (e). For perpendicular 
polarization, the field D is smooth. We then write 



(E) = (e- 1 D) = (e- 1 )(D) 



(46) 



and e±_ = (e -1 ) -1 . 

Similar considerations can be applied to the 2D prob- 
lem of Sec. IIIIDI For waves polarized along the z-axis, 
the field E is smooth, which results in e zz = (e), in agree- 
ment with ([44)) . 

One can also consider a more general smooth field of 
the form S = pE + (1 — p)T> = [p + (1 — p)e]E, where 
p is a mixing parameter. Here we consider the 3D case 
and assume that S is smooth for any polarization state. 
Application of the smooth field principle results in the 
following equalities: 



(E) = (S)(l/[p+(l-p)e]) 
(D) = (S) (e/\p+(l-p)e]) 



(47a) 
(47b) 



from which we find the effective permittivity to be 



(e/[e + p/(l-p)]) 

(l/[6 + p/(l-p)]> 



(48) 



Eq. ([48]) is, in fact, the Maxwell- Garnett formula. Al- 
though this form is rarely used, the Maxwell- Garnett ef- 
fective permittivity can be written as 



<e/(e + 2e 6 )) 
(V(c + 26 6 )) 



(49) 



We see that (ggj and (g9j) coincide if p = 2e b /(l + 2e&). 

Thus, the Maxwell- Garnett EMT assumes that the 
field S = [(e + 2e 6 )/(l + 2e 6 )]E is smooth. Since the 
mixing parameter p depends on the permittivity of the 
host medium, Eq. ([49]) is not invariant with respect to 
the substitution e a O e& and p O 1 — p. The homog- 
enization formula ([33]) derived in this paper, however, 
is fully symmetric. Note that Bruggeman's EMT is also 
symmetric but can not be easily written in terms of aver- 
ages. Therefore, it is not clear which form of the smooth 
field Bruggeman's approximation assumes. In general, 
the smooth field does not need to be a linear functional 
of E and D. 



IV. REFLECTION AND REFRACTION AT A 
HALF-SPACE BOUNDARY 

An infinite lattice is a mathematical abstraction. All 
experimental media are bounded, and the physical ef- 
fects which occur at the boundary are often important. 
For instance, as mentioned above, it is not possible to 
determine simultaneously and uniquely the effective per- 
mittivity and permeability of a medium from the bulk 
dispersion relation alone. 
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The problem of reflection and refraction of a wave at a 
flat interface is considered in this section. The goals are 
three- fold. First, we will derive the limit in which the 
correct expression for the Fresnel reflection coefficient is 
obtained. This will turn out to be the same limit as was 
used in Sec. IIII Al Second, we will show that the correct 
expression for the reflection coefficients results only if we 
take £ = 1 in ([32]) . from which it follows that Jl = 1. 
Third, we will provide additional mathematical justifica- 
tion for the results of Sec. IIII A[ Indeed, the derivations 
of that section contain one dubious step. Namely, the 
applicability of the Poisson summation formula ([H]) can 
be questioned because the variable q is complex. Strictly 
speaking, the series in the left-hand side of ([T4j) diverges 
for an infinite lattice. The problem can be fixed, in princi- 
ple, by considering real- valued q's and then analytically- 
continuing the summation result to the whole complex 
plane. In this section, no such complication will arise 
since all series in question are convergent. 



x 























n 













































z 



-oo < n x ,n < oo 
1 < n <oo 



FIG. 2: (color online) Sketch of the geometry considered: re- 
flection and refraction at a half-space boundary. 



A. General setup 

The geometry considered in this section is sketched in 
Fig. [21 The medium occupies the right half-space and 
the left half-space has the background permittivity q,. It 
would be more appropriate to consider the case when the 
left half-space is vacuum and the right half-space is a two- 
component mixture, so that there are three different com- 
ponents in the problem. This, however, requires the use 
of the half-space Green's tensor Q - a step that is not 
conceptually difficult, yet mathematically involved. Here 
we restrict consideration to only two components. This 
includes the cases when the host medium is vacuum and 
also when the incident beam is first refracted from vac- 
uum into a homogeneous medium of permittivity ^ 1 
(at a planar interface that is located at z = z\ <C — ft 
and is not considered explicitly) and then into a het- 
erogeneous medium which is a mixture of a- and 6-type 
components. 

Physically, the z coordinate of the effective medium 
boundary can be stated only approximately, within an 
interval of width ~ ft. It will prove mathematically con- 
venient to place the boundary on the plane z = 0, and 
the centers of the left-most cells on the plane z = ft, 
as shown in Fig. [2] In the EMT developed below, the 
half-space z > is assumed to be filled with an effective 
medium. 

A wave can not propagate in a semi-infinite medium 
without an external source. Therefore, we must solve 
the integral equation ([2]) with a nonzero incident field E$ 
which we will take to be a plane wave. We will also find 
that, under appropriate conditions, a uniquely-defined 
reflected plane wave E r exists in the region z < 0. The 
incident and the reflected waves are given by 

E^(r) = Ai exp(lq • r) , — oo < z < oo , (50a) 
E r (r) = A r exp(k r • r) , — oo < z < . (50b) 



Note that the incident wave is defined in the whole space 
but Eq. (j2j) is only defined for r G ^tot- The wave num- 
bers of the incident and the reflected waves can be written 
as 

k^ = kj_ -|- zki z , k r = k_L zkiz • (^-0 

Henceforth, the subscript "_L" will be used to denote pro- 
jections of vectors onto the x?/-plane. Note that k^ -z = 
and 

k\ = k 2 r = k 2 ± + kt = k 2 b = k 2 e b = (^) 2 e b . (52) 

It is important to note that the vector k^ is purely real. 
A complex- valued k^ would imply a wave that is evanes- 
cent in a direction parallel to the interface. This would 
necessitate the presence of additional interfaces; such a 
possibility is not considered here. The vector k^ is real- 
valued even if the host medium is absorbing. Indeed, 
we should keep in mind that the incident wave enters 
the host medium from vacuum and that the tangential 
component of the wave vector is conserved at any planar 
interface, even if one of the media is absorbing. How- 
ever, the z-projection of k^ does not need to be real. In a 
transparent host (e& > 0), the incident wave is evanescent 
and k iz is purely imaginary if k± > k^; in an absorbing 
host, ki Z is, generally, complex. 

Note that the reflected wave ()5Qbj) does not enter 
Eq. ([2j) because it is identically zero in fitot- The reflected 
wave is computed a posteriori once the polarization field 
P is found. Then the amplitudes A r and A$ can be used 
to determine the reflection coefficient. 

To solve Eq. ([2|) in the presence of the incident field, 
we decompose P as 

P = P B + Ps , (53) 

where P# is the Bloch wave of the form (jSJ) and P^ is 
an additional wave that originates due to the presence of 
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the surface. We seek the condition under which 



E EO (r)= I G(r,r')P B (r' 

J r2tr>t 



)d 6 r 



E B (r)+Eext(r) + E s (r) , (54) 



where in fL 



E B (r) = ^P B (r) , 
E ext (r) = -Ei(r) , 

If (|55 1) -(jS1) l) hold, then Eq. (0) becomes 



(55a) 
(55b) 



Ps(r) 



3x 

47T 



E 5 (r)+/ G(r,r')P s (r')rfV 



r G fitot ■ (56) 



Note that Eq. ([56]) contains only quantities which are 
associated with the surface wave. 

Eq. ([54]) is the mathematical formulation of the Ewald- 
Oseen extinction theorem and we will refer to Eeo as to 
the Ewald-Oseen field. We will see that one can deter- 
mine the reflection coefficient from the conditions ([55]) . 
We will also see that the surface wave is exponentially 
localized near the interface and does not contribute to 
either reflection or transmission coefficients if 



(ki 



^f>kt 



Vgx^O 



(57) 



Inequality ([57]) is weaker than what is required for ho- 
mogenization. It is merely the condition that there is 
no Bragg diffraction in the medium; if ([57]) is violated, 
the conventional reflection and transmission coefficients 
can not be defined. If, however, ([57]) holds, we do not 
need to solve Eq. ([56]) explicitly; it suffices to know that 
the surface wave does not contribute to any measurement 
performed sufficiently far from the interface. 



B. Evaluation of the Ewald-Oseen field 

To compute the Ewald-Oseen field, we proceed along 
the lines of Sec. lIII Al to arrive at the following expression: 

EEo(r) = l/^ p) I d3i?F ( R) 

x exp [ip • (r - R)] ^ exp [i(q - p) • r m ] . (58) 

m 

So far, no restrictions on r have been placed. In partic- 
ular, r can be either in the right or left half-space. How- 
ever, when we later substitute the result of integration 
into Eqs. ([55]) . r will be restricted to f^tot- 

The sum over m in ([58]) can be evaluated as follows. 



First, we expand the summation as 
^exp [z(q-p) -r m ] = 

m 

oo 

^ exp [i(q x - p x )hm x + i(q y - p y )hm y } 



m x ,m y = — oo 
oo 

x ^ ex.p[i(q z - p z )hm z ] . 

m z =l 



(59) 



From symmetry considerations, we know that = k^. 
This property is a manifestation of momentum conser- 
vation and will be confirmed below by considering the 
conditions ([55]) . Since, as discussed above, kx is purely 
real, q x and q y are also real. Therefore, we can com- 
pute the sums over m x and m y using the Poisson sum 
formula. Further, the half-range sum over m z converges 
absolutely because the transmitted wave decays into the 
medium and, correspondingly, Im^ > 0. We, therefore, 
have 

^exp [z(q-p) -r m ] 

m 



/(Pz)$>(px-qx-g±) , (60) 



where 



/(^) = Yl exp ( q * Pz "> hrriz \ 

m z =l 

1 



ex.p[i(p z - q z )h] - 1 
2tt ( 27r _1 



h ^ p z -q z - 9z 



(61a) 



(61b) 



Here the well-known Laurant expansion of the function 
l/[exp(zz) — 1] has been used. The equality (|61b|) is an 
important observation. It will allow us to evaluate the 
Ewald-Oseen field. 

We now proceed by substituting ([60]) into ([58]) . which 
yields 

E EO (r) = ^E /_~ l^f(P*)K(qi_ + gx + zp z ) 

x / d 3 i?F(R)exp[i(q ± +g± +zp z ) ■ (r-R)] . 
Jo. 

(62) 

The integral over p z can be computed by contour inte- 
gration since all the poles and residues of the integrand 
are known. The positions of the poles in the complex 
Pz-plane are shown in Fig. [3l The poles at p z = q z + g z 
are the singularities of the function f(p z ). Since q z has a 
positive imaginary part and all g^s are real- valued, these 
poles lie in the upper half-plane. The remaining poles are 
the singularities of K(q± + + zp z ), which is viewed 
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lmp z 



p z = A 2 -(q±+g±) 2 «fe± 



P Z = & - 



Re/? z 



FIG. 3: (color online) Poles of the integrand of Eq. (|62[) in 
the complex p 2 -plane. 



here as a function of p z . From the definition (j6j) , we find 
that these singularities are located at p z = ±V S± , where 



k 2 b ~ (qx + gx) 2 • 



(63) 



These poles can be considered separately for g± = 
and gx 7^ 0. The t wo poles corresponding to gx = 
are p z = ±Vo = ±\/^6 — ^he P°l es w ^ n g± 7^ 
have large (either positive or negative) imaginary parts 
if hkb,hq_\_ <C 1, in which case they can be written, ap- 
proximately, as V s± ~ ig±- 

Note that in the case of infinite lattices, the singulari- 
ties of K(p) do not contribute to Fourier integrals of the 
type ([T2]) because the corresponding residues are identi- 
cally zero [these singularities fall in between the peaks of 
the delta-function fence given by the right-hand side of 
dll|) ]. We will compute the contributions of the different 
families of poles to the integral ([62]) separately. If the 
vector of position r is inside one of the inclusions, the 
integration contour must be closed in the upper half of 
the complex p^-plane. Correspondingly, only the poles 
with positive imaginary parts contribute to the integral 
([62]) in this case. The Ewald-Oseen field can also be 
computed in the left half-space. If the point of obser- 
vation r is further away from the interface than ft/2, so 
that the inequality z r < — ft/2 holds, the integration 
contour must be closed in the lower half of the complex 
p^-plane. In what follows, it will be shown that the poles 
at p z = q z +g z yield the Bloch-wave field E^(r), the pole 
at p z = Vq yields the extinction field E ext (r), the pole 
at p z = — Vo yields the reflected wave, and, finally, the 
poles p z ~ i^g^ with gx 7^ yield the fast-decaying 
surface wave. 



1. Block wave 

We start by computing the Bloch-wave contribution 
to the Ewald-Oseen field, E#(r). We place the point 
of observation in fitot? use the expression (|61b|) for f(p z ) 
and evaluate the contributions of the poles p z = q z + g z to 
the integral ([62]) . This results in the following expression: 

Eb W = S ex P [* (q + g) • r] i^(q + g) 



/ F(R)exp[-z(q + g) -H}d 3 R , r G ft t 



(64) 



Here we have used the equalities gx + zg z = g and 
J2 S± ^Z Sz = Zlg- Now, if F(R) is expanded accord- 
ing to ([T6]h and if the expansion coefficients F g satisfy 
(fT7|). then the field given by Eq. ([64)) satisfies E B (r) = 
(47r/3x)Ps(r) for r G ft to t, where P# is of the form (|gj). 
Thus, (|55a|) is satisfied if the Bloch wave of the polar- 
ization Pb is the same as one would find by solving the 
eigenproblem ([T7]) for an infinite medium. This justifies 
the use of the Poisson summation formula in Sec. IIII At 
Eq. ([T7]) applies to general photonic crystals that are 
not necessarily describable by EMPs. As was discussed 
in Sec. IIII A[ homogenization is obtained by taking the 
limit ft — >• 0. This limit must be computed separately for 
the equations with g = and g ^ 0, which results in ([T9]) . 
This system of equations defines an eigenproblem for the 
Bloch wave vector q, while the polarization vector Fo 
is obtained as an eigenvector of ([22]) . The higher-order 
expansion coefficients F g are uniquely determined by Fo 
but Fo itself is defined by ([T9]) only up to a multiplicative 
factor. Next, we will show that this factor is fixed by the 
condition ([55b)) . 



2. Extinction wave 

We now compute the contribution of the pole located 
at p z = Vo- The function f(p z ) is analytic in the vicin- 
ity of Vo] therefore, we can use the expression ([6 lap for 
f(p z )- Since Eqs. (|55b|) should hold only for r G ^tot> we 
close the integration contour in the upper half-plane. A 
straightforward calculation yields 



Eext(r) = 



47ri exp [i (qx + zVo) • r] 



ft 2 exp [i (Vo - q z ) ft] - 1 

h b 



kl - (qx + ZVo) <8> (q_L + z^o) 



2Vo 

x / d?KF{K) exp [—i (qx + zVo) • R] , (65) 

r G fitot • 

We seek the condition under which E ext (r) = — E^(r) 
for r G fitot? where E^(r) is given by ([50ap . It immedi- 
ately transpires that the above equality can hold only if 
qx = kj_. The continuity of the tangential components 
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of all wave vectors, including the incident wave vector Iq, 
the reflected wave vector k r and the Bloch wave vector of 
the transmitted wave q follows from the discrete transla- 
tional symmetry of the problem. We now find from ([5T]) 
that Vo = ki Z and + zVo = . With the use of these 
equalities and the notation 



F(k)= / F(R)exp(-ik-R)d*R 



(66) 



we can simplify Eq. ([65]) as 

4m exp (iki • r) 



Eext(r) 



h 2 exp[i(k iz -q z )h] 



2ki 



-P(ki) , r G O t , 



The extinction condition then takes the form 
2m k 2 -ki® ki F(k^) 



(67) 



(68) 



h 2 exp [z (fc^ - q z ) h] - 1 fe^ 

So far, no approximations have been made. The homog- 
enization limit is obtained by observing that 

lim exp [i (±k iz - q z ) h] = 1 + i(±k iz - q z )h , (69a) 

h—^Q 

lim F(ki) = lim F(k r ) = V(l + S)F . (69b) 

Once the above limiting expressions are used, the extinc- 
tion condition becomes of the form 

This equation couples the amplitude of the incident field, 
A^, and the amplitude of the Bloch polarization wave, 
Fo- The vector Fo must simultaneously satisfy the fol- 
lowing two conditions: (i) be an eigenvector of the tensor 
in the square brackets in Eq. ([22]) and (ii) satisfy ([70]) . 
These two conditions determine both the direction and 
the length of F . 



4- Surface wave 

Finally, let us evaluate the contribution of the poles 
Pz = P gl with gx 7^ 0. For r G ^tot? we nave 



k b k s± g k g 



^F(kg^) , rG^tot (73) 



where 



(74) 



If the condition ([57]) holds, the quantities T^g^ have 
nonzero imaginary parts even if the host is transparent. 
Therefore, the surface wave decays exponentially away 
from the interface. In the homogenization limit, the ex- 
ponential decay is fast. Indeed, in the limit h —± 0, we 
have (for g ± ^ 0): V s± -> ig±, k s± -> g ± + izg±, 



-l/g±h. With these limits taken into ac- 



count, the surface wave takes the following form: 

2m ^ k 2 - (g ± + zz#_l) <8> (g ± + zz#_l) 



/i 3 



x exp [(ig ± - zg±) • r] F (g ± + zg±) , (75) 

r G fitot • 

It can be seen that E5 decays exponentially on the scale 
of h. So does the wave of polarization P#, as both fields 
are related by the integral equation ([56]) . 

Solving Eq. ([56]) numerically can be a very difficult 
task. Fortunately, doing so is not necessary if one is only 
concerned with far-field measurements. 



3. Reflected wave 

Consider now the case when the point of observation r 
in the left half-space. As discussed above, we will place 
r at least h/2 away from the interface. This will allow us 
to close the integration contour in ([62]) in the lower half 
of the complex p^-plane. The reflected wave is obtained 
by computing the input of the pole p z = — Vo- We find 
that the electric field of the reflected wave is of the form 
(|50bp where the amplitude A r is given by 

A = 2m k 2 -k r ® k r F(k r ) 
h 2 exp [-i (ki Z + q z )h]-l k iz 

This expression contains no approximations. In the ho- 
mogenization limit, we use the limiting expressions ([69]) 
and obtain 

A — 4ra (1+s > F °- <72) 



C. Reflection coefficient 

We will now utilize the results of the previous subsec- 
tion to compute the reflection coefficients for the half- 
space. We will use the assumption of Sec. IIII Bl namely, 
that the crystallographic and optical axes of the medium 
coincide so that the tensor E is diagonal in the labora- 
tory frame. Apart from other simplifications, media of 
this type are non-chiral and do not rotate the polariza- 
tion of the transmitted and reflected waves. This prop- 
erty holds even beyond the homogenization limit, since 
it is a straightforward consequence of the elementary cell 
symmetries, and it will enable us to consider the s- and 
p-polarizations separately. 

In this subsection, we will explicitly use the reference 
frame shown in Fig. [2j That is, we will assume that the 
plane of incidence is the xz-plane and that the projection 
of the wave vectors k^, k r and q onto the interface is 
k_L k x ~x.. 
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1. S -polarization 

In the case of s-polarization, the incident and reflected 
waves are polarized perpendicularly to the plane of inci- 
dence. Consequently, we have A^, A r oc y, and the exact 
reflection coefficient is given by 



A r • y _ F(k r ) • y exp [i(k iz - g z )h] - 1 
Ai • y P(ki) • y exp [-i(k iz + q z )h] - 1 



. (76) 



To derive the second equality, we have used the expres- 
sions ([68]) and ([7T]) for the amplitudes A^ and A r . This 
is an exact expression that retains its physical meaning 
as long as ([57]) holds. In the homogenization limit, we 
use the expressions ([69]) to obtain 



Here q z is given by 



/x» iy^\ Is* £ 



(77) 



(78) 



which follows from the dispersion relation ([37]) . in which 
we must take q x = k x . The square root branch in ([78]) is 
determined by the condition lm(q z ) > 0. 

The expressions ([77]) and ([78]) should be compared to 
the corresponding Fresnel coefficient rp and the disper- 
sion relation for a homogeneous medium characterized by 
the permittivity and permeability tensors e and ft: 



ki. 



r F 



(79) 



The wave number q z in an effective medium satisfies the 
dispersion relation 



q z 



k 2 t 



yyf^xx 



h2 ^xx 
x • 



(80) 



As was discussed in Sec. lIIIB"T| we must impose the con- 
dition ([32]) on the EMPs e and ft in order to obtain the 
same laws of dispersion in the composite and in the con- 
tinuous medium models. In particular, this condition 
guarantees that the quantities q z given by Eqs. ([78]) and 
([80]) are equal for all values of k x . But if this is so, the 
only way the two expression ([77]) and ([79]) can yield the 
same reflection coefficient is if we set £ = 1 in ([32]) , which 
corresponds to ft = 1. 

We note that to reach the above conclusion, it is 
sufficient to consider the reflection coefficient for s- 
polarization only. We will show next that the same con- 
clusion can be reached by considering p-polarization only 
and that the homogenization results obtained in these 
two cases are consistent. 



2. P-polarization 

In the case of p-polarization, the reflection coefficient 
can be conveniently defined by using the ratio of tangen- 
tial components of the magnetic field for the reflected and 



incident waves. The magnetic field amplitudes of these 
waves are given by 



(81) 



As could be anticipated, the amplitudes B^ r are aligned 
with the y- axis. We can now use the expressions ([68]) and 
([7T]) for the amplitudes A^ r to find the exact reflection 
coefficient: 



>r • y 



B 7 



•y 



k r x F(k r ) 


' y exp [i(k iz - 


q z )h] - 1 


k» x F(ki) 


. ^ exp [-i(k iz - 


\-qz)h\-\ 



(82) 



In the homogenization limit, this expression is simplified 
by using (|6U|) . which leads to 



[kr x (1 + E)F ] • y k iz - q z 
' [ki x (1 + S)F ] • y k iz + q z 



(83) 



As shown in Appendix [Bj Eq. ([83]) can be further simpli- 
fied to read 



kj z /e h - q z /r] x 
h z /e h + q z /rj x 



(84) 



In ([83 ]) , ([M]) , q z satisfies the dispersion relation for the 
p-polarized wave, ([38]) . With the substitution q x = k x , 
the latter reads 



k 2 r] x - k 2 



(85) 



As in the case of s-polarization, the branch of the square 
root is determined by the condition lm(q z ) > 0. 

We wish to compare the expressions ([84]) and ([85]) to 
the analogous expressions in a continuous medium with 
the EMPs e and ft. The Fresnel reflection coefficient for 
a p-polarized incident wave is given by 



r F 



ki Z /e h - q z /e a 



ki Z /e b + q z /e xx 
and the dispersion relation in the effective medium is 



(86) 



k 2 e x fi y k 2 - 



(87) 



As in the case of s-polarization, the condition ([32]) with 
an arbitrary parameter £ guarantees that the two expres- 
sions ([85]) and ([87]) yield the same wave number q z for 
all values of k x . However, the expressions ([84]) and ([86]) 
yield the same reflection coefficient only if we set £ = 1 
in (p2j). 

This completes the proof that the correct choice of 
the parameter £ in ([32]) is £ = 1 and, correspondingly, 
the correct homogenization result is ft = 1. A similar 
proof has been given by us for a one-dimensional layered 
medium in Ref. [43| for both s- and p-polarizations. 
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V. COMPARISON OF POINT-DIPOLE AND 
CONTINUOUS-MEDIUM MODELS 

The model of point-like polarizable particles arranged 
on a three-dimensional infinite lattice possesses an intu- 
itive physical appeal. Historically, many authors have 
used this model and, although an exhaustive review is 
outside of the scope of this paper, Refs. (334361 Eol EM |46| 
can be mentioned. Unfortunately, the model is haunted 
by divergences. In this section, we will discuss the na- 
ture and origins of these divergences and some of the 
commonly-used methods for their regularization. We will 
also attempt, to the degree it is possible, to establish a 
correspondence between the model of point dipoles and 
the model of a continuous two-component medium, which 
is the subject of this paper. 

Most previous works on electromagnetic waves in 
point-dipole lattices assume that the background 
medium is vacuum. For compatibility of results and sim- 
plicity of notations, we will also make this assumption (in 
this section only) and set = 1 -MO, kb = k = uo/c + iO. 

The model of point dipoles considers an array of point- 
like particles which have well-defined locations, but no 
shape or size. Instead of the latter two physical charac- 
teristics, the electric dipole polarizability a is used. In 
some generalizations of the model, the magnetic dipole 
polarizability is also included. The basic idea of this ap- 
proach is that the electromagnetic response of a particle 
is completely characterized by its polarizability. 

If only the electric polarizability is retained, one ar- 
rives, in lieu of the integral equation (|2]), at the set of 
algebraic equations 



—d ri 

a 



Ei(r n )- 



G(r n , r^dn 



Here d n is the electric dipole moment of the n-th particle. 
Now two important points should be made. First, the 
summation on the right-hand side of ([55]) is restricted 
only to the indices m which are not equal to n. This 
reflects the idea that the electric field at the site of the nth 
dipole is a superposition of the incident wave E^(r n ) and 
the waves scattered by all other dipoles. Second, energy 
conservation requires that [47l-[49| Im(l/a) < —2k 3 /3. If 
the equality holds, the particles are non-absorbing. It is 
convenient to decompose the inverse polarizability as 



1 

a 



1 



2k 3 



(89) 



where «ll is the "Lorenz-Lorentz" quasistatic polariz- 
ability and —i2k 3 /3 is the first non- vanishing radiative 
correction to the imaginary part of 1/a. Radiative cor- 
rections to the real part of 1/a also exist and are, in fact, 
of a lower order in /c, but it is the correction to the imag- 
inary part which is physically important and should be 
retained even in the limit kh —> 0. We will see momen- 
tarily that the two seemingly unrelated facts mentioned 
above are mathematically connected. 



We now consider an infinite lattice, set the incident 
field to zero and seek the solution to ([55]) in the form 
d n = dexp(zq • r n ). This results in the eigenproblem 



-d = 5(q)d 

a 



where 



(90) 



0] (91) 



is the dipole sum. Using the Fourier representation (j5j), 
we rewrite (I9T1) as 



s(q) 



47T 

Y 



/ 



(2tt) 2 



K(p) V] exp [i(p - q) • (r„ 



(92) 

The first complication encountered in the above is that 
the summation on the right-hand side of ([92]) is incom- 
plete. We can easily fix this problem by adding and sub- 
tracting unity to the series, which leads to 



5(q) 



47T 

T 



d 3 p 
(2^)3 



K(p) 



(93) 



where we have used the Poisson summation formula ([H 
Still, both terms on the right-hand side of ([93]) are diver- 
gent. We will deal with the integral first. To this end, 
we utilize the expression for K(p) given in (j6j) and notice 
that the angular integral of the term p 2 — 3p <g> p is zero 
in three dimensions. Therefore, we have 



d 3 p 
(2iry 



K(p) = 4tt 



p 2 dp 2k 2 



(2tt) 3 p 2 



k 2 



(94) 



This is still a divergent integral. We can regularize ([94] 
by writing 



= lim I in [ 
^ I Jo 



p 2 dp 2k 2 
(2tt) 3 p 2 - k 2 



exp[-(Ap) 2 



(95) 



The above limit indeed exists and is equal to i/c 3 /27r, 
assuming that Imk > (which is true if we take k = 
u/c + zO). Upon substitution of this result into ([93]) . we 
find that 



47T 

T 



2k 3 



(96) 



We now use the decomposition ([59]) and notice that the 
above term is canceled by a similar term on the left-hand 
side of ([90]) . Taking into account this cancellation, ([90]) 
becomes 



g 



(97) 



The mathematical tricks used so far are not very objec- 
tionable. The result ([96]) is a reflection of the fact that 



lim 



4ttA 3 



/ 

J\r' 



G(r,r')dV 



|r'-r|<A 



2k 3 



(98) 
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Here we have assumed that the particle is spherically 
symmetric. The use of a different integration volume 
in ([98]) . or of a different regularization function in ([95]) . 
would certainly yield a different result. Fortunately, if 
kh <C 1, only the real part of / is affected by the choice 
of the regularization function in ([95]) while the imaginary 
part is relatively stable. If Re/ is unimportant, e.g., if 
it is small compared to the sum of real parts all other 
contributions in ([93]) . then ([97]) is a good approximation, 
regardless of the true shape of the particles. 

However, the divergence of the series in the right- 
hand side of ([97]) is truly problematic. One can attempt 
to regularize this divergence by the same mathematical 
trick that was used above. However, the result of such 
a manipulation would indeed depend on the regulariza- 
tion function in a nontrivial way. One can conclude that 
knowledge of the particle polarizability is, in fact, insuf- 
ficient for solving the problem at hand. The shape of the 
particles is also important and can not be disregarded. 

Another way to look at this is the following. The po- 
larizability a defines the response of a particle to an ex- 
ternal electric field which is almost uniform over the par- 
ticle volume. However, in an infinite three-dimensional 
lattice, the electric field is not uniform over the parti- 
cle volume, no matter how small the particle is. This is 
because the lattice Green's function W(R, R/) given by 
([12]) experiences an integrable divergence when R = R/. 
However, in the point-dipole model, we are attempting 
to evaluate this function exactly at R = R' = 0, which 
is not mathematically reasonable. 

It appears that the only feasible approach to regular- 
ize the summation in ([97]) is to endow the particles with 
a finite volume, as was done, for example, in Ref. [34j . 
This would naturally lead to a modification of ([97]) in 
which the right-hand side is multiplied by a decaying 
function /(g), ensuring convergence. Unfortunately, the 
exact form of /(g) strongly depends on the particle 
shape and size. If the regularization is carried out in 
a mathematically-consistent way, one would end up with 
a set of equations that are identical to the equations ob- 
tained here, for the model of a continuous two-component 
medium. 

Evidently, within the point-dipole model, one wishes 
to avoid introducing the particle shape and size. Then 
the only conceivable approach to regularization is simply 
to truncate the series in ([97]) . by leaving only the g = 
term in the summation, which leads to the eigenproblem 



1 47T 

— d=—K(q)d. 



(99) 



Regularization of this type is, in fact, appropriate for 
small spherical particles. If one also uses the quasistatic 
polarizability of a sphere of radius a, namely, 

1 



3 c a 



(100) 



then ([99]) becomes equivalent to the Clausius-Mossotti 
relation and the EMT that follows from it is the standard 
Maxwell- Garnett approximation. 



One may be tempted to forget about the limits of appli- 
cability of Eq. ([99]) . In other words, once ([99]) is derived, 
it is technically possible to use it with any polarizability 
c^LL- The latter can be obtained independently, i.e., by 
solving the Laplace equation for a single isolated parti- 
cle of arbitrary shape. Unfortunately, this approach is 
mathematically inconsistent. Eq. ([99]) was derived from 
([97]) by applying a regularization method which is only 
appropriate for small spheres. Applying ([99]) to particles 
of nonspherical shape is likely to result in errors. 

In summary, the model of point dipoles is capable of re- 
producing the standard Maxwell- Garnett mixing rule for 
small spheres. Radiative corrections to this result can 
also be derived (35) . However, in three dimensions, the 
model breaks down and can not be used when a substan- 
tial deviation from the Maxwell- Garnett approximation 
is expected, i.e., for particles whose volume fraction is not 
small or whose shape is different from a sphere. In other 
words, the model does not provide a mathematically con- 
sistent way of computing the self-energy E which appears 
in equations ([23]) or ([33]) and is, therefore, usable only in 
the physical situations when £ can be neglected. Never- 
theless, we note that in systems of lower dimensionality 
(e.g., in nanoparticle chains), the point-dipole model is 
useful and can provide significant physical insights. 



VI. CONTINUED-FRACTION EXPANSION OF 
THE SELF-ENERGY AND THE MEAN-FIELD 
APPROXIMATION 



A. Abstract notation 

In this section, we will find it convenient to rewrite 
Eqs. (]19b|) and ([2T]) in Dirac notation. First, we note 
that, in order to recover all components of the tensor £, 
one must solve ([19b)) for three different right-hand sides: 
Fo = x, Fo = y and Fo = z. To this end, we introduce a 
triplet of infinite-dimensional vectors \ap), operators Q, 
M, W, and vectors \bp) (/? = x,y,z) according to 



(ag\ap) = M(g)5 a £ 



(ag\Q\a'g f ) 
(ag\M\a'g f ) 



' (1 - ^9cx9cx') 

'M(g-g') , 



W = QM , 
\bp) = Q\ap) . 



(101a) 
(101b) 
(101c) 
(lOld) 
(lOle) 



Note that Q is diagonal in the index g, M is diagonal in 
the index a, but the product of the two, W = QM, is 
not diagonal. We must also keep in mind that the index 
g in the above equations is not allowed to take the zero 
value. We further define the vectors \Fp) as the solutions 
to 



±-w)m = \bp). 



(102) 
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The above is equivalent to the set (|19b|) . The tensor 
elements of E are defined by 

= (a a \F p ) = (a*\ (— - w) \b ) 
V PX J 



(103) 



(o a |( QM) Q\a p ) 



It can be seen that E is computed as the resolvent of the 
operator W = QM and plays the role of the self-energy, 
which accounts for interactions between the inclusions. 



B. Mean- field approximation 

The mean-field approximation is often misunderstood. 
In particular, it is unrelated to Maxwell- Gar nett theory. 
Rather, it allows one to replace certain operators by ap- 
propriately chosen scalar multiples of the identity. The 
approximation reproduces the exact zeroth and first mo- 
ments of the resolvent and serves as the first-order ap- 
proximation in its continued-fraction expansion. Here 
the a ppr oximation is explained following Berry and Per- 
cival [50|. 

Let us seek the solution to Eq. ([102)) in the form 
\Fp) = A 1 6/3), where A is a scalar to be determined. Upon 
substitution of this ansatz into (|102)h we obtain the equa- 
tion 



1 1 

PX A 



\h) = w\h) • 



(104) 



Because |&#) is, generally, not an eigenvector of W, there 
is no such value of A for which Eq. (|104p would hold. 
The best we can hope for is that a projection of this 
equation onto a given vector would hold for some A. Since 
we are interested not in the whole vector \Fp) but in 
its projection onto \a a ), it seems reasonable to project 
Eq. (fTMl) onto the latter. This yields 



A 



PX 



1 - PX 



(Oa\W\h) 

(a>a\bp) 



(105) 



and the corresponding mean-field approximation for the 
self-energy is 



^a(3 — 



PX(aa\bp) 



pxMQWp) 



(a a \W\bp) 1 (a a \QMQ\a/3) ' 
1 - PX—, — rr^— 1 - PX~ 



(a a \Q\a(3) 



(106) 



As was mentioned in Sec. IIII C[ the matrix element 

(a a \Q\a p ) = [M(-g)Q(g)M(g)]^ (107) 

is identically zero for inclusions with cubic symmetry (in 
three-dimensional composites) so that Eq. (|106|) yields in 



this case zero and is not useful. If (a a \Q\a^) is zero, a 
non- vanishing mean-field approximation can be obtained 
by "shifting" the solution according to \Fp) = px\bp) + 
\Fp). The self-energy is then given by Eq,^ = (a a \Fp) 
where \F*) satisfies 



PX 



W)\F'}=pxW\bp) . 



(108) 



The mean-field approximation for the "shifted" equation 
flTUgD is 



{px) 2 (a a \QMQ\ap) 



1 - PX 



(a a \(QM) 2 Q\a p ) 
(a a \QMQ\a ) 



(109) 



C. Continued- fraction expansion of the self-energy 

Continued- fraction expansions (CFEs) are very useful 
in physics [5l|, [52[. The mathematical underpinning of 
all CFEs is the theory of the correspondence between the 
formal Laurent series of meromorphic functions and cer- 
tain continued fractions [53) . There exists a deep math- 
ematical relation between CFEs and the problem of mo- 
ments, that is, the problem of finding a distribution from 
the knowledge of its moments. 

CFEs can be derived in different ways. Hay dock (5l| 
has employed the Lanczos recursion to transform a cer- 
tain Hamiltonian to tridiagonal form. A diagonal ele- 
ment of the inverse of a tridiagonal matrix can be writ- 
ten as a J-fraction (a continued fraction of Jacobi type). 
In Ref. [51], this procedure was applied to a Hermitian 
operator to compute a diagonal matrix element of the re- 
solvent. In this paper, the operator W in ([102)) or ([103)) 
is not symmetric or Hermitian and we are interested in 
off-diagonal elements of the resolvent. Therefore, the nu- 
merical procedure used by Haydock is not directly appli- 
cable. Perhaps, it can be generalized to become appli- 
cable to the problem at hand; we have not explored this 
possibility. Instead, we will derive a CFE for the right- 
hand side of Eq. ([103)) from the following theorem which 
does not require any symmetry properties of the opera- 
tors involved, yields a CFE for arbitrary off-diagonal el- 
ements, and, to the best of our knowledge, has not been 
reported in the literature. The resultant expansion will 
be an S-fraction (a continued fraction of Stieltjes type). 
Note that an S-fraction can always be transformed into 
a J-fraction by the so-called equivalence transformation. 

Theorem 1 Let W be a linear operator acting on the 
Hilbert space T~L and Z be a complex number. Suppose 
that |^> e H. If (i) + and (ii) (Z - W)' 1 

exists, then 



(<p\(z-wy 



{(p\{Z- WT)~ l W\ip) 



(110) 
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where 



an) 



a) Circle 



b) Square 



The proof is given in Appendix^ Note that hll(J\) has a 
finite limit when Z — » 0. 

The factor {4>\{Z — WT)~ 1 W\iIj} in the denominator of 
([TTQj) can be written as {<j>\(Z - Wx)' 1 ^), where Wi = 
WT and = W\ip). The formula (fTTUj) can now be 
applied to Wi) -1 and so on iteratively. After 

some manipulation, this yields the following expansion: 



{<j>\{Z - 



^2 



(112) 



1 



^3 



Note the interlacing factors of Z and 1. The coefficients 
Kj (j = 1,2,...) are obtained from a three-point recur- 

= and 



sion. Namely, starting from = 0, 
Ki = (4>\^)j we compute for j = 1,2,... 



iv J -+i) = w(iv J ->-« j |^_i>) , 



Kj+i - 



(4>\4>j) 



(113) 



To obtain a CFE of the right-hand side of Eq. (|103|h we 
identify Z 1/px, W QAL \<j>) \a a ) and 
\bp)=Q\a p ). 

With the above substitutions taken into account, it 
transpires that the coefficients tv j are determined only 
by the geometry of the composite. Once a set of Kj have 
been found for a given geometry, the EMPs can be easily 
computed for any material parameters of the composite 
constituents. This is a characteristic feature of a spectral 
theory and the CFE (|112p is, in fact, a spectral represen- 
tation of the self-energy E. 

Finally, we note that, in the case of three-dimensional 
composites with cubic symmetry, the first condition of 
the Theorem does not hold when the theorem is applied 
directly to (|103p . In this case, one can build a CFE 
staring from the "shifted" equation (j!08|) . 



VII. NUMERICAL SIMULATIONS 

A. General setup 

Numerical simulations have been performed for a two- 
dimensional composite. The composite is periodic in the 
xy-plane while the inclusions form infinitely-long fibers 
which are oriented parallel to the z-axis and can have 
different cross sections. The case when the electric field 
is parallel to the fibers is not considered here, since this 
polarization results in a simple arithmetic average of the 
type dSj) . However, when the electric field is polarized in 
the xy-plane, the homogenization problem is nontrivial 
and can be numerically challenging. We will consider 
inclusions with circular and square cross sections, as is 





FIG. 4: (color online) Two types of elementary cells used in 
numerical simulations. 



illustrated in Fig.HJ The functions M(g) for these shapes 
are given in Appendix lAl 

It is assumed that the host medium is vacuum and the 
inclusions are metallic and characterized by a frequency- 
dependent Drude permittivity of the form 



= 1 



3c4 



uj(uj + ij) ' 



66 = 1 • 



(114) 



In Eq. (|114p . uop = u p /\^3 is the Frohlich frequency, 
uj p is the plasma frequency, and 7 is the Drude relax- 
ation constant. We will compute the effective permit- 
tivity of the composite e as a function of frequency for 
0.1 < uj/ujf < 2 and for the fixed ratio j/ujf — 0.1. It is 
assumed that, for all frequencies used in the simulations, 
the basic condition for the validity of a standard EMT, 
kbh,qh <C 1, is satisfied. 

Numerical simulations will be performed by truncating 
the infinite set of equations (|19b|) so that the vectors g 
fill the box 



2nL/h < g x ,g y < 2nL/h 



(115) 



where L is an integer. The total number of g- vectors 
which satisfy the above inequality is (2L + l) 2 and the 
total number of algebraic equations to be solved is N = 
2[(2L + l) 2 — 1], where we have accounted for the fact 
that the vector g = is excluded in the set of equations 
(|19b|) . It can be seen that N — > SL 2 when L —> 00. In 
the simulations, we will use integer powers of 2 for L, 
up to L = 2 8 = 256. The latter case corresponds to 
N = 526, 366 equations. 

The truncated set of equations (|19b|) can be solved by 
any direct numerical method. The computational com- 
plexity of direct methods is 0(N 3 ) and the solution must 
be obtained anew for every frequency used (we sample 
the frequency at 200 equidistant points in the interval 
0.1 < u/uf < 2). This is time-consuming but possi- 
ble for L < 64. For larger values of L, direct methods 
become impractical. We will use, therefore, the CFE 
of Sec. IVI CI The computational complexity of this ex- 
pansion is 0(jmax^ 2 ) 7 where j max is the order of trun- 
cation of the continued fraction. More specifically, the 
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continued fraction is truncated by assuming that kj = 
for j > j max , so that only the first j max coefficients are 
used in Eq. (j!12|) . For the problem at hand, j max « 50 
will prove sufficient. Other iterative methods, such as 
the conjugate gradient method, also have computational 
complexity 0(j max 7V 2 ), j max being the number of iter- 
ations. However, the computationally-intensive part of 
the conjugate-gradient solver (when applied to Eq. (|19b|) ) 
must be repeated for every value of lj, while the coeffi- 
cients Kj in (|112p need to be computed only once for a 
given geometry. 

The inclusions shown in Fig. [4] have cubic symme- 
try. As was discussed in Sec. IIIIB1 the self-energy E 
is reduced in this case to a scalar. As a result, the ef- 
fective medium is isotropic in the x?/-plane. Of course, 
anisotropy can still be revealed if the polarization vector 
has a component along the z-axis. In the simulations re- 
ported below, we have computed E by solving Eqs. (|19b|) 
and using the definition (|2T|) . The effective permittivity 
for transversely-polarized waves was then computed by 
using Eq. (j33j). 



B. Convergence and stability 

The convergence of the CFE (|112|) with the trunca- 
tion order of the continued fraction, j max , is illustrated 
in Fig. [5] Here the real and imaginary parts of the effec- 
tive permittivity are plotted as functions of frequency. It 
can be seen that the convergence is very fast for circular 
inclusions and somewhat slower for square inclusions. In 
all cases, jmax — 50 is sufficient for convergence. 

The three-point recurrence relation (|113|) is numeri- 
cally unstable for large values of j. This is illustrated in 
Fig. [6] Shown in this figure are the coefficients k,j ob- 
tained on two different computers for the geometry de- 
scribed in the figure caption. The same code and input 
data were used in both cases. The coefficients from the 
two sets coincide for j < 50 with high precision. However, 
differences start to appear at j ~ 50 and, at j ~ 100, the 
coefficients are unreliable. The instability occurs when 
an iteration step in ([113)) asks for a relatively small dif- 
ference of two large numbers and the numerical precision 
of the floating-point arithmetic is exceeded. 

The instability illustrated in Fig. [6] appears to be trou- 
blesome but is, in fact, of little concern. This is illus- 
trated in Fig. [7J which displays the effective permittivity 
computed by the CFE (jl 12j) for various truncation or- 
ders jmax, and the same quantity computed by solving 
Eqs. (|19b|) directly. One of the sets of kj's displayed in 
Fig. [6] has been used for computing the data points for 
panels (a,b) of Fig. [71 Despite the instability, the curves 
with j max = 50 and jmax = 100 are indistinguishable 
and very close to the data points obtain by direct inver- 



sion of (|19bp . Thus, the unreliable coefficients Kj do not 
influence the final result. This is one of the nice proper- 
ties of all CFEs: a numerical instability does not result 
in numerical imprecision. It is true that increasing the 
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FIG. 5: (color online) Convergence of the CFE (fTT2l) with the 
truncation order j max for circular (a,b) and quadratic (c,d) 
inclusions with the same volume density p — 0.16. The set of 
equations (|19b|) has been truncated using L = 64. In panels 
(a,b), the curves with j max = 30,40, 50 are indistinguishable. 



truncation order beyond j max = 50 is not useful, but it 
is not harmful either. This point and some related issues 
are discussed in more detail in Sec. IVIIll below. 

Having established the convergence properties of the 
CFE, we next consider convergence with the size of the 
box, L (up to now, all plots have been computed for L = 
64). In Figs. \S}9\ e is plotted as functions of frequency 
for various values of the density, p, and the box size, L. 
Also shown in these figures are the results obtained from 
the generalized Maxwell- Garnett formula 
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FIG. 6: (color online) Absolute values of the coefficients Kj 
computed for circular inclusions with p — 0.16 and L — 64 on 
two different computers (CI and C2). The same FORTRAN 
code and input data have been used in both cases. 



2p 



3 e b + u(e a - e h ) 

P (fa ~ €b) 

3e b + v(e a - e b ) 



(116) 



which applies to ellipsoids, v being the appropriate de- 
polarization factor. In the case of three-dimensional 
spheres, v = 1/3 and Eq. (|116p coincides with Eq. ([33]) 
in which the self-energy E is set to zero. In the case 
of infinite circular cylinders, the depolarization factor, 
which corresponds to the orthogonal electric polarization, 
is v = 1/2. 

Several conclusions can be drawn from Figs. [8] and 
First, convergence is obtained for boxes of reasonable 
size. In all cases shown, L = 256 yields very accurate 
results, and in some cases L = 64 is sufficient. However, 
it is important to note that we have verified the conver- 
gence by doubling the size of the box. Determination of 
convergence by using linearly sampled values of L, (say, 
L = 10,11,12...) can be misleading. This is a typical 
situation when boundary- value problems are solved nu- 
merically. Convergence must be established by at least 
doubling the size of the mesh used. 

Second, it can be seen that convergence is faster for 
p = 0.32 than for p = 0.16. Although the electromag- 
netic interaction is stronger in the second case, the faster 
convergence is to be expected. Indeed, the size of the box 
should be selected so that the sum rules (|A2|) are satis- 
fied with some reasonable precision, and that is achieved 
at smaller values of L for larger values of p. Even faster 
convergence is obtain for p = 64 (data not shown). How- 
ever, at the percolation threshold (p = 7r/4 & 0.79 for 
circular inclusions), the convergence is relatively slow. 

Third, the generalized Maxwell- Gar nett formula (j!16p 
with v = 1/2 yields a reasonable result for circular in- 
clusions with p = 0.16. Even better agreement has been 
obtained for p = 0.08 and p = 0.04 (data not shown). 



3 
2 
1 


-1 

4 

3 
2 
1 



1 1 


1 


Re(e) f\ 


- 


(a) Circles 




1 1 


^ 1 



Im(e) 



(b) Circles 




1 

Panels a-d: 

J max — 10 
jmax = 50, 100 

Direct 



1 1 

Re(e) 


1 


— (c) Squares IV, 






FIG. 7: (color online) Effective permittivity of circular (a,b) 
and square (b,c) inclusions computed using the CFE (|112[) 
with the truncation orders j max = 10, 50, 100, and by direct 
inversion of Eqs. (|19b|) . In all cases, p — 0.16 and L = 64. 
The data points for j max = 50 and jmax = 100 are visually 
indistinguishable and, therefore, represented with the same 



However, as the size of circular inclusions increases, the 
Maxwell- Garnett approximation becomes less accurate. 
For square inclusion, the approximation is inaccurate 
even for very small values of p. In all cases, the electro- 
magnetic interaction tends to shift the absorption peaks 
from the Maxwell-Garnett's prediction towards the lower 
frequencies. At p = 0.32, the effect is already quite pro- 
nounced. 
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FIG. 8: (color online) Convergence of the effective permit- 
tivity e with the size of the box, L, for circular (a,b) and 
square (c,d) inclusions with p = 0.16. The curves labeled as 
v — 1/2 and v — 1/3 have been obtained from the general- 
ized Maxwell- Garnett mixing formula (|116|) for the values of 
v indicated. 



C. Comparison of inclusions of various size 

We finally compare the effective permittivity for circu- 
lar and square inclusions of different sizes. The results are 



displayed in Figs. 1101111 In the case of circular inclusions, 
there exists a pronounced spectral peak which shifts to- 
wards lower frequencies when p is increased. However, 
once the inclusions touch (this happens at p = 7r/4 w 
0.79, the single resonance is destroyed and a broad ab- 
sorption band develops. The lower- frequency behavior 
of e is in this case metallic, since the percolating sample 
is characterized by a nonzero static conductivity. This 
result can not be obtained within the Maxwell- Garnet 
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FIG. 9: (color online) Same as in Fig. [8] but for p — 0.32 and 
different values of L, as indicated. 



approximation, or the Bruggemann approximation, even 
at a qualitative level. 

The square inclusions do not touch for p < 1. Corre- 
spondingly, the low-frequency behavior of e is not metal- 
lic even for large filling fractions, e.g., for p = 0.85. In- 
terestingly, at relatively small values of p, the absorp- 
tion spectrum forms a band with one main resonance 
and many minor resonances which are shifted towards 
the shorter waves. However, as p increases, the minor 
resonances become less pronounced. At p = 0.85, the 
spectrum is dominated by a single Lorentzian-type res- 
onance. In the case of circular inclusions, the picture is 
somewhat different. A single Lorentzian resonance ex- 
ists at small values of p and additional minor resonances 
develop as p increases. These additional resonances are 
clearly visible in the p = 0.64 curve shown in the left 
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FIG. 10: (color online) Effective permittivity for circular in- 
clusions of different volume densities. The p ~ 0.79 case cor- 
responds to the percolation threshold (touching circles). 



column of Fig. [TUl 



VIII. DISCUSSION 

A few points that deserve additional discussion are ad- 
dressed in this section, in no particular order. 



A. Conditions of applicability 

The EMT derived in this paper describes a composite 
medium accurately if qh,kh <C 1. There are no addi- 
tional conditions. In particular, there is no requirement 
that the permittivity (or conductivity) of any constituent 
of the composite be bounded. However, if a metallic in- 
clusion has very small losses (very high conductivity), 
then the effective permittivity computed according to the 
formulas of this paper can have one or more sharp spec- 
tral peaks. These peaks are caused by electromagnetic 
resonances in the inclusions (which we have not disre- 
garded by any means) and can be seen in Figs. ISlfTTl In 
the spectral regions where these resonances take place, it 
is possible that q ^> k. This effect is known as the res- 
onance wavelength shortening. Conceivably, the Bloch 
wave number q can become so large due to this effect 



FIG. 11: (color online) Same as in Fig. \W\ but for square 
inclusions. 



that the condition qh <C 1 would break. In this case, our 
theory is inapplicable. 

The above consideration can be construed as a justifi- 
cation for development of extended EMTs, e.g., by taking 
a limit in which the conductivity of metal inclusions goes 
to infinity first [23|, [24], 127I. or by using other trajectories 
in the parameter space [26|. However, two important 
caveats exist. First, in many known applications, EMPs 
of the order of unity are required, e.g., e ~ ft « — 1 is 
required for operation of a superlens. In this case, of 
course, q ~ k, there is no resonant wavelength shorten- 
ing, and our theory applies. The second caveat is that, 
even if metal inclusions have very high conductivity, the 
imaginary part of the obtained effective permittivity is 
not small close to a resonance. This can be clearly seen 
in Figs. [MTU Therefore, there is not much hope to obtain 
a resonant effect without having, simultaneously, strong 
absorption in the medium. This observation is in agree- 
ment with Stockman [54] , although we do not pursue here 
a rigorous mathematical consideration of this point. 

Finally, in the case when qh is not actually small com- 
pared to unity and our theory does not apply, it appears 
from considering the exact reflection coefficients ([T6]h([82]) 
that any EMPs that can be introduced in any theory 
would depend on the angle of incidence. More generally, 
the EMPs would depend on the type of illumination. We 
conclude that the medium is simply not electromagneti- 
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cally homogeneous in this case. 



B. The case of small losses 

Another problem associated with high conductivity of 
metallic inclusions is numerical stability and convergence. 

The simulations of Sec. I VIII have been performed for a 
relatively large loss parameter, j/ujf = 0.1. If this num- 
ber is substantially reduced, the convergence with the 
truncation order of the continued fraction, j ma x, is ex- 
pected to become slower. A general rule of thumb is that 
the truncation order should not be less than the number 
of clearly discernible peaks in the function lme(uj) (the 
absorption spectrum). This is because the CFE trun- 
cated at the order j max captures correctly the first j max 
moments of the above function. At sufficiently large val- 
ues of j, the three-point recursion ([113)) becomes numer- 
ically unstable, as is illustrated in Fig. [6j If the required 
value of jmax is larger than the value of j at which the 
onset of numerical instability occurs, then the CFE will 
not yield an accurate numerical result. 

The situation outlined above is common for all iter- 
ative methods. For example, the convergence of the 
conjugate-gradient method becomes extremely slow for 
small ratios of j/ljf] at some point, the recurrence re- 
lations used in the conjugate-gradient iterations also be- 
come numerically unstable. One can hope to improve 
stability by noting that the n-th order tail of the CFE 
(j!12p . that is, the expression 



is also an expansion of a certain resolvent, and the in- 
stability occurs because the parameter e (defined in the 
proof of Theorem 1, Appendix IC]) becomes numerically 
small. This can be fixed by "shifting" the operator A as 
described in Sec. IVIB1 In this way, a nested set of CFEs 
can be obtained, where each CFE is numerically stable, 
as well as the whole expression. 



C. Consideration of chirality and polarization 
conversion 

Although the general formalism of this paper allows 
one to take chiral media into consideration, all deriva- 
tions which were brought to a logical conclusion have 
been carried out for the non-chiral case. This has pro- 
vided a mathematical simplification, yet left untouched a 
wealth of interesting physical phenomena which are asso- 
ciated with chirality. This shortcoming will be addressed 
by us in the future. 

Even if the medium is non-chiral, it can exhibit the 
effect of polarization conversion [55|, which has been re- 
cently predicted and experimentally observed in deeply- 
subwavelength nanostructures in Ref. [13|. In Sec. IIV C\ 



we have made an assumption that the plane of incidence 
coincides with one of the crystallographic planes of the 
medium. In this case, the s- and p-polarized waves are 
independent and polarization conversion does not occur. 
However, the homogenization result obtained in this pa- 
per is more general and, in particular, it is applicable to 
any direction of incidence. If the plane of incidence does 
not coincide with any crystallographic plane, the geome- 
try of the problem becomes similar to that considered in 
Ref. [13| and polarization conversion can occur. In other 
words, the reflected and transmitted (in the case of a fi- 
nite slab) waves due to a purely s- or p-polarized incident 
wave can have both s- and p-polarized components and, 
at least theoretically, it is possible to design a medium 
with the conversion coefficient close to unity. 



D. 3D vs 2D simulations 

So far, we have performed simulations only for 2D me- 
dia. One can argue that in the 3D case the size of the 
algebraic problem would become so large as to render 
the method unusable. Of course, three-dimensional elec- 
tromagnetic problems are always challenging. However, 
there is reason for optimism. Namely, the formula for 
the effective permittivity ([33]) uses the three-dimensional 
Maxwell- Garnett approximation as the point of depar- 
ture. In other words, a nonzero value of E provides a 
correction to the three-dimensional Maxwell- Garnett for- 
mula. This happens to be true even for two-dimensional 
media. However, the three-dimensional Maxwell- Garnet 
formula is inaccurate in the 2D case even for very thin 
cylinders, as is clearly illustrated in Figs. I8|91 In the nu- 
merical simulations of Sec. IVIII (for circular inclusions), 
a lot of effort was spent to compute accurately the self- 
energy E whose effect was, essentially, to transform the 
Maxwell- Garnett from a 3D to a 2D form. 

In the case of small three-dimensional inclusions, one 
can expect a much faster convergence with L. For ex- 
ample, if the inclusions are small spheres, an accurate 
result is obtained by starting with E = 0. As the spheres 
increase in size, the Maxwell- Garnett approximation be- 
comes less accurate and a nonzero value of E must be 
used. However, as we have seen in the numerical simu- 
lations, the required values of L are, in fact, smaller for 
larger sizes of the inclusions. 

Mathematically, the above considerations are related 
to an interesting fact which was mentioned in Sec. IVII 
Namely, the matrix element (a a \Q\ap) is identically zero 
for three-dimensional cells with cubic symmetry. Conse- 
quently, the mean- field approximation and the continued- 
fraction expansion must be derived for the "shifted" 
equation (|108j) . As a result, the mean- field formula (|109|) 
contains an overall factor of (px) 2 while in the 2D simu- 
lations of Sec. I VIII this factor was equal to px- 
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IX. SUMMARY 

We can draw the following conclusions: 

1. A medium constructed from nonmagnetic compo- 
nents is also nonmagnetic in the limit h — >• 0. This 
result is in line with arguments put forth in (28^ . 
the simulations in [29j and the more formal math- 
ematical theory of [25j. 

2. The model of point-like polarizable particles is ill- 
suited for homogenization of three-dimensional pe- 
riodic composites due to inherent divergences. The 
point-dipole approximation can be still a useful the- 
oretical tool for studying systems in lower dimen- 
sions. 

3. In agreement with the previous conclusion, we have 
found numerically that the EMPs are sensitive to 
the shape of inclusions even if the volume frac- 
tion is small. Thus, circular and square inclusions 
in Figs. I7|8I have very different spectra of EMPs, 
even though the volume fraction of the inclusions 
is p = 0.16. When the volume fraction becomes 
larger, the differences between the circular and the 
square shapes are dramatic. Thus, it is shown in 
Figs. II 01 lTl that the percolation phenomenon occurs 



for the circular inclusions at the volume fraction 
p = 7r/4 « 0.79, when the inclusions touch. The 
composite in this case is conducting. The compos- 
ite consisting of square inclusions of the volume fill 
fraction (which do not touch) is still a dielectric. 

4. We believe that the goal of homogenization theory 
is to describe a given physical composite. There- 
fore, rather than studying different limits, which 
correspond to different trajectories in the param- 
eter space, it is important to delineate regions of 
the parameter space and to determine, to which 
one of these regions the particular composite be- 
longs. Along similar lines, we note that a satis- 
factory theory of homogenization requires error es- 
timates. That is, it is critical to understand how 
the error in the homogenization limit depends upon 
contrast. We plan to investigate this question in fu- 
ture work. 
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VM(-g)M(g) = - , 
g p 

VM(g-g / )M(g / ) = -M(g) . 

f J n 



(A2b) 
(A2c) 



These equations hold for inclusions of arbitrary shape. 
Now define a complimentary function iV(g) by 



^/ expHg-R)**. 



(A3) 



Appendix A: Mathematical properties of M(g) and 
some special cases 



Here C denotes the unit cell and C\fl is the region com- 
plimentary to the inclusion. It can be seen that iV(g) has 
all the properties of M(g) with the substitution p — >> 1—p. 
Additionally, the functions M(g) and N(g) are related 

by 



P M(g) + (1 - p)N( S ) = 5 e0 . (A4) 
From this, we obtain the low and high-density limits: 



lim 7V(g) lim M (g) S s0 . (A5) 
p— >-0 p— >i 



Of course, the high-density limit is unreachable for most 
regular shapes (with the exception of cubes). For exam- 
ple, in the case of spheres, the maximum allowed value 
of p is 7r/6. 

Some special cases of M(g) are given below. For an 
inclusion in the shape of either a 3D sphere or a 2D circle 
of radius a < h/2, 



From the definition ([T8|h it follows that 



M 3D (g) = 



M 2D (g) 



3[sin(ga) — gacos(ga)] 



(9a) 3 



2J1M 



ga 



(A6a) 
(A6b) 



M(0) = 1 , M(-g)=M*(g) . 



(Al) 



For the case of inclusions whose center of symmetry coin- 
cides with the center of the unit cell, we have M(— g) = 
M(g) and, therefore, M(g) is real. If the center of sym- 
metry is displaced by a vector a, the function M(g) is 
transformed according to M(g) — >• exp(— ia. • g)M(g). 

By applying the Poisson summation formula, we can 
derive the following sum rules: 



where J\{x) is the cylindrical Bessel function of the first 
kind. For a parallelepiped or rectangle centered at the 
origin with all faces parallel to the crystallographic planes 
and sides of length 2a Xl 2a y and 2a z , 



M 3D (g) = 



M 2D (g) 



sin(g x a x ) sm(g y a y ) sin(g z a z ) 

9x&x Qy^y 9z&z 
sm(g x a x ) sm(g y a y ) 



9y a y 



(A7a) 
(A7b) 
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Appendix B: Details of some calculations pertaining 
to the case of p-polarization 

To simplify notations, we will denote (in this Appendix 
only) 



1 + £ = S , px = ^, 



so that 



1 + 2KS a a 
1 — K>S a a 



(Bl) 



(B2) 



We start by deriving Eq. ([39]) . To this end, we write the 
wave vector of the p-polarized wave as q = q x it+q z z (note 
that q y =0) and seek a nontrivial solution to Eq. ([22]) . 
Multiplying ([22]) by the non-zero factor q 2 — k 2 and using 
(|6]), we obtain the following equation: 

(q 2 - k 2 )F - n(2k 2 + q 2 )SF + 3^q (q • SF ) = . (B3) 

We now account for the fact that the tensors E and S = 
1 + E are diagonal in the laboratory frame and write 



and 



(SF ) a = S aa F 0a , a = x,y,z 



q • SFq — q x S xx Fo x + q z S zz Fo z . 



(B4) 



(B5) 



Using this result, and projecting Eq. ([B3[) onto the y- 
axis, we immediately obtain Fq v = 0. The two remaining 
Cartesian components of Fo satisfy a system of two linear 
equations, which are obtainable by projecting (|B3|) onto 
the x- and z-axes. These two equation are not linearly 
independent, provided that the dispersion relation ([38]) 
holds [otherwise, the only solution to ([B3[) is trivial]. It 
is, therefore, sufficient to consider one of these equations, 
say, by projecting ([B3[) onto the x-axis. The resultant 
equation is 

AF 0x + BF 0z = , (B6) 

where 

A = (1 - kS xx ) q 2 + 3kS xx & - (1 + 2kS xx ) k 2 b , (B7a) 



B = 3nS zz q x q z . 



(B7b) 



We now simplify the expression (|B7aj) for the coefficient 
A. Specifically, we substitute into this expression q 2 — 
q 2 +q 2 and k 2 = e^k 2 = e\ ) (q 2 /r] x J rq 2 /r] z ), where we have 
used the dispersion relation ([38]) . This yields 



A = (1 - kS xx ) (q 2 z + q 2 x ) + 3f^S xx q 2 x 
-e b (l + 2KS xx ) + ^ 



(B8) 



We now use (|B2ft to write out the quantities r] x and rj z 
in (|B8|) in terms of S xx and S zz . It can be seen that the 
terms proportional to q 2 cancel, and we obtain 



a q o 1 + 2nS xx 2 

i± — oKjD zz — — — q 
1 + 2kS^ k 



(B9) 



(BIO) 



We use this result and the expression (|B7b[) for B to 
compute 

Fox^ B^ 1 + 2/c5 gg 
F 02 A 1 + 2^5^^ 

Returning to the original notations ([Bl|h we obtain ([39]) . 

Next, we show how to derive Eq. ([84]) from ([83]) . 
Eq. ([83]) contains the factor 

[k r x(l + E)F ]-y _ [k r x5Fo]-y 
" [ki x (1 + E)F ] • y [k 4 x 5F ] • y ' 1 ; 

which we will now evaluate. To compute the vector prod- 
ucts, we note that Iq = xfc^ + zfc^, k r = itk x — zfc^ and 
SFq — &S xx Fq x + zS zz Fq z . From this, we find 



kxSzzFo z + ki z S xx F( 



Ox 



k>x&zzFo z ki Z S xx Fo x 



(B12) 



Next, we use the ratio Fq x /Fq z given by (|B10|) . account 
for the conservation of the wave vector projection onto 
the interface, that is, q x = k X: and re- write (]B12|) as 



R 



k 2 S zz (1 + 2kjS xx ) ki Z q z S xx (1 + 2k,S zz ) 
k 2 S zz (1 + 2kjS xx ) + k{ Z q z S xx (1 + 2kS zz ) 



. (B13) 



To proceed, we need to exclude the variable k 2 from 

(|B13|) . Using the dispersion relations ([38]) and ([52]) for 

the refracted and the incident waves (in the geometry 
considered, q 2 = k\ — k 2 ), we write 



/7 2 P 1 1 

- h — — ft — — K b — — \K X -h AC i; 



6fe 



6fe 



(B14) 



Solving (|B14|) for /c 2 , we obtain 



l/»7* " 1/C6 
1 + 2kS zz ( 1 — K>S XX ^ t 2 



1 + 2/sS x 



(B15) 



where we have used <\B2\i to obtain the second expres- 
sion from the first. We now substitute the result given 
in ([BT5]) into flBT3j) . The factors of 1 + 2kS zz in the 
numerator and the denominator cancel, and we obtain 



R 



(1 hzS xx )q 2 (1 + 2KjS xx )k 2 z 3tvS xx ki z q z 



(1 — KjS xx )q 2 — (1 + 2tvS xx )k 2 z + 3tvS xx ki Z q z 

(B16) 

At the next step, we divide the numerator and the de- 
nominator in (jB16|) by the factor 1 + 2kS X x and, account- 
ing for the identity 



1 + 2kS x 



1 1 

£6 Vx 



(BIT) 



obtain 



n 2 k 2 

Vz ^iz 



R = 



1 



kizQz 



Vx Cb \Cb Vx 

Qz k 2 z (I 1 . ? 

— - — + kizQz 

Vx £b \£b Vx, 



(B18) 
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The expressions in the numerator and denominator can 
now be factorized, and we arrive at the final result 



(Qz + hz) 



R 



kiz Qz 



(Qz - hz) ( — + ^ 



(B19) 



Substitution of this expression into (f83j) immediately re- 
sults in ([84]) . 



Appendix C: Proof of Theorem 1 
C.l. An equivalence transformation 



To derive the equality (jllO|) . we first introduce some 
notation. Let 



a = 



Zl - {^(s + BP^B^) ' 



(C5) 



C.2. A useful identity 

Below, we will frequently use the following identity: 

{(t)\B\^) =Z<r-e. (C6) 
The above equation is easily derived by noting that 

(<p\B\^) = {4>\{z-w)- l w\^) 

= (4>\(z - w)~\w - zm + z(4>\(z - wy 1 ^) 

= -e + Za . (C7a) 



£ = 


, 


(Cla) 


P = 




(Clb) 


R{Z-A) = 


(Z-A)- 1 , 


(Clc) 


B 


R(Z;W)W , 


(Cld) 


a = 


(<t>\R(Z;W)\1>) ■ 


(Cle) 



Here R{Z\ A) is the resolvent of the linear operator A 
and Z is a complex number. In the new notation, the 
operator T defined in (|1 1 lj) takes the form 



T = 1 - -P 

€ 



and Eq. ([110p is rewritten as 
1 



Z 1--(0|^(Z;WT)WW) 



(C2) 



(C3) 



Note that, by the first hypothesis of the Theorem, e ^ 0. 

We now write the following chain of equalities in which 
the second hypothesis of the Theorem, namely, that 
R(Z; W) exists, has been used: 

r(z-wt) = (z-wry 1 

= (z- W +-WP^J 

'r- 1 (Z;W) + jWP^J 

R~\Z\ W) ^1 + \r{?\ W)WP > j 
= e[e + R(Z; W)WP]~ 1 R(Z; W) . (C4) 



Using the last equality in (|C4p and the notation (|Cld|h 
we rewrite (|C3p identically as 



C.3. The main derivation 

To proceed, we need to express the operator (e + 
BP) -1 , which appears in the right-hand side of (|C5[) . 
in a more tractable form. To this end, consider the equa- 
tion 



(s + BP)\x) = \b) , 



(C8) 



where \x) is viewed as the unknown and \b) ^ is an 
otherwise arbitrary element of the same Hilbert space. 
Using the definition of P (|Clb[) . we transform (|C8|) to 

e\x)+B\<l>)(<l>\x) = \b) , (C9) 
project the result onto and find that 



BY 



(CIO) 



We now use the previously-derived identity (|C6|) in the 
right-hand side of (|C10p to obtain 



Za 



(Cll) 



Upon substitution of (jCllj) into ()C9|) . we find the solution 
to jnHl) or ([C9]), namely, 



e v -z^J |6) • (C12) 

Since the vector \b) in (|C8|) is arbitrary, we conclude that 



( e + flP)-i = I(l-*!gM) . (C13) 
This equality can be verified directly by substitution. 
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C.4. Putting everything together 



We can now put everything together and obtain (|C5|) . 



From HCm . we have (<j>\(e + BP^B^) = ^' |y/ = 1 - — . (C15) 

m + «P)- = | (i - « M = 1 , (Ci4) Upon s _„ ot this result int0 the righ , hand 5ide 



of (|C5|) . we find that the latter is, indeed, an identity, 
where we have, again, used (|C6|) . Now, we can write and so are (|C3j) and (jllOp . 



